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Background  and  Scope: 

Mid-frequency  (1-10  kHz)  shallow  water  acoustics  is  of  important  Naval  interest,  but  had  been 
an  area  where  concentrated  basic  research  was  relatively  lacking  before  the  year  2000.  This 
might  have  been  a  consequence  that  while  long-range  propagation  (low  frequency,  100’s  of  Hz) 
dominated  ASW  work,  direct-path  (high  frequency  >10  kHz)  found  applications  in  MCM, 
leaving  mid-frequency  research  being  covered  under  either  program,  but  without  sufficient 
emphasis.  As  an  outgrowth  of  APL-UW's  high-frequency  tradition,  we  started  planning  mid¬ 
frequency,  shallow  water  research,  approaching  the  subject  by  emphasizing  contemporaneous 
environmental  and  acoustics  studies  in  order  to  quantitatively  understand  dominant  physical 
mechanisms.  The  major  components  of  the  research  are  transmission  loss,  bottom  and  sub¬ 
bottom  impact  on  propagation  and  scattering,  surface  and  water  column  influence  on 
propagation,  culminating  in  a  full  understanding  of  the  physics  of  shallow  water  mid-frequency 
reverberation.  Under  this  project,  we  conducted  a  series  of  theoretical  and  numerical  studies, 
supplemented  by  existing  data,  to  plan  for  a  comprehensive  shallow  water  reverberation 
experiment,  which  resulted  in  the  TREX13  field  project. 

Modeling  shallow  water  reverberation  is  a  problem  that  consists  of  two-way  propagation 
(including  multiple  forward  scatter)  and  a  single  backscatter.  In  order  to  understand  the 
reverberation  problem  at  the  basic  research  level,  both  propagation  and  scattering  physics  need  to 
be  properly  addressed.  Some  aspects  of  reverberation  are  better  treated  stochastically  — 
roughness  scattering  from  the  bottom  and  surface,  for  example.  Other  aspects  can  be  more 
successfully  treated  deterministically,  such  as  waveguide  propagation  and  target  scattering.  Yet 
there  are  aspects  where  experience  is  limited,  and  the  approach  used  to  model  them  would  likely 
be  determined  based  on  further  empirical  investigations.  Scattering  from  fish  schools  and 
shipwrecks  along  with  biological  aggregates  around  them  falls  into  this  category.  Clutter  is  a 
major  issue  and  often  unexplainable  by  known  objects  in  the  waveguide.  We  proposed  a  new 
hypothesis  for  clutter:  Some  of  the  clutter  is  due  to  the  combination  of  (1)  forward  scattering  of 
propagating  sound  from  low-grazing  angle  to  high-grazing  angle,  and  (2)  backscattering  of  the 
high-grazing  angle  energy,  creating  a  target-like  clutter  because  high-grazing  angle  backscatter  is 
much  greater  than  low-grazing  angle  backscatter. 

Although  there  had  been  reverberation  measurements  at  various  frequencies,  there  had  not  been  a 
true  6. 1  level  reverberation  experiment  where  the  environment  has  been  sufficiently  measured  to 
support  full  modeling  of  the  data.  It  is  only  after  all  components  that  contribute  to  reverberation 
are  well  measured  and  modeled  that  a  true  understanding  of  the  reverberation  problem  can  be 
achieved.  During  this  funding  period,  we  concentrated  on  making  preparations  for  the  TREX13 
experiment,  which  included  both  theoretical  development  and  data  analysis. 
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Results 

The  main  result  of  is  work  manifested  in  the  successful  field  experiment,  TREX13.  While 
TREX13  had  a  broad  design  scope  and  covered  almost  all  aspects  of  shallow  water  acoustics,  a 
central  theme  based  on  a  SONAR  equation  as  given  in  Fig.  1  summaries  the  philosophy  and 
general  approach.  The  SONAR  equation  appropriate  for  the  project  is: 

RL  =  SL  -  2x  TL  +  ISS 

where  RL  is  reverberation  level,  SL  the  source  level,  2xTL  the  two-way  transmission  loss,  and 
ISS  the  scattering  strength  integrated  over  the  scattering  patch  for  given  sonar  beam.  Two  unique 
features  of  TREX13  are  that:  1)  all  components  of  the  SONAR  equation  are  designed  to  be 
individually  measured  in  the  same  frequency  band  over  the  same  environment,  2)  an  extensive 
environmental  measurements  at  the  appropriate  temporal  and  spatial  resolutions  would  be  made 
such  that  basic  research  questions  concerning  predictability  and  uncertainty  of  shallow  water 
reverberation  can  be  quantitatively  addressed.  The  approach  to  data  analysis  can  be  summarized 
into  the  following  steps: 

1.  Based  only  on  acoustic  measurements,  assess  to  what  degree  the  measured  reverberation, 
transmission,  and  scattering  quantities  satisfy  the  SONAR  equation?  This  first  step 
establishes  a  complete  data  set  that  enables  detailed  follow-up  analysis.  It  also  bounds  the 
predictability  and  uncertainty  of  reverberation. 

2.  Incorporating  environmental  data,  assess  the  predictability  and  uncertainty  of  the 
individual  terms  in  the  SONAR  equation.  Identify  key  environmental  parameters  that 
contribute  to  the  variability. 

3.  Review  available  reverberation  predictive  models,  and  if  necessary,  develop  new  models 
to  incorporate  environmental  knowledge  in  order  to  improve  model  accuracy  and/or 
speed. 

4.  With  both  acoustics  and  environments  measured,  divide  model  environmental  parameters 
into  categories,  e.  g.,  those  that  can  be  inferred  from  acoustic  data  and  those  where 
databases  are  necessary. 

5.  Given  prediction  requirements  and  uncertainty  tolerance,  provide  a  set  of  key 
environmental  parameters  necessary  as  input  to  models. 
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Why  this  equation? 


Related  Websites 


TREX13 

Reverberation  with  Supporting  Environment 


RL  =  SL  -  2x  TL  +  ISS 


Figure  1.  Schematic  of  TREX13  based  on  TREX13  website  under  development.  The  400  kHz 
backscatter  strength  from  the  multibeam  survey  (de  Moustier)  is  shown  as  the  gray  swath  along 
the  Main  Reverberation  Track  where  focused  environmental  characterization  was  conducted. 
Main  acoustic  assets  are  shown  here:  The  reverberation  sources  and  receiving  arrays  are 
deployed  from  the  R/V  Sharp  and  are  labeled  “Sharp,  Reverb”.  The  three  vertical  line  arrays 
fielded  by  Scripps  Institution  of  Oceanography  are  shown  as  triangles  with  labels  given  distances 
to  the  reverberation  source. 

Some  of  the  specific  results  documented  in  the  literature  and  are  summarized  as  follows: 

1.  Shallow  water  reverberation  modeling.  A  key  question  for  evaluating  reverberation 
models  is  speed  vs.  fidelity.  Without  fidelity,  a  model  would  not  be  able  to  make 
appropriate  prediction;  without  speed,  the  model  could  not  be  effective  in  applications. 
During  this  funding  phase,  we  developed  a  Monte  Carlo  reverberation  model  [2]  based  on 
PE  and  perturbation  theory  which  allows  incorporating  specific  environmental 
knowledge  while  still  maintaining  reasonable  computational  speed.  The  parabolic 
equation  method  is  used  to  handle  the  two-way  propagation,  and  first  order  perturbation 
theory  is  used  to  handle  the  backscatter.  Because  the  calculation  time  is  independent  of 
the  number  of  realizations,  this  method  is  much  faster  numerically  than  any  models 
available.  Another  advantage  of  this  method  is  that  it  can  easily  handle  complications 
such  as  internal  waves  and  swells. 

2.  Clutter  is  often  considered  to  come  from  large  scatterers.  We  hypothesized  that 
reverberation  clutter  can  often  result  from  a  combination  of  forward  scatter  and 
backscatter.  As  an  example,  we  investigated  the  case  where  non-linear  internal  waves  can 
change  sound  propagation  path  to  high  grazing  angles,  and  those  high  angle  sound  would 
in  turn  be  backscattered,  resulting  in  a  ‘ghost’  clutter’  [1],  This  line  of  investigation  will 
continue  in  the  analysis  phase  of  the  TREX13  data. 
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3.  Transmission  Loss  is  one  of  the  major  issues  in  shallow  water  acoustics  because  it  is  a 
quantity  against  which  detection  predictions  are  made.  A  couple  of  papers  [5,  6]  were 
devoted  to  this  topic  where  influences  of  oceanographic  variability  are  considered. 

4.  How  to  measure  through  either  direct  method  or  inversion  techniques  is  a  key  to  a 
successful  for  a  6.1  field  project  because  understanding  the  environmental  impact  on 
acoustics  is  the  foundation  for  any  quantitative  analysis.  A  couple  of  papers  [3,7]  address 
the  environmental  issues,  one  concerns  the  inversion  of  geo-acoustic  parameters  from 
chirp  sonar  data,  the  other  documents  progress  of  a  piece  hardware  which  can  take  direct 
measurements  on  sediment  sound  speed  and  dispersion.  Both  techniques  have  been 
further  developed  during  TREX13. 

Now  that  the  TREX13  field  work  has  been  accomplished,  the  next  phase  is  to  analyze  the 
collected  data  and  document  the  results.  The  anticipated  areas  of  advances  are  shallow  water 
transmission  loss  and  reverberation  clutters,  both  are  areas  having  basic  research  and  application 
needs. 
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Abstract:  Clutter  is  related  to  false  alarms  for  active  sonar.  It  is 
demonstrated  that,  in  shallow  water,  target-like  clutter  in  reverberation 
signals  can  be  caused  by  nonlinear  internal  waves.  A  nonlinear  internal 
wave  is  modeled  using  measured  stratification  on  the  New  Jersey 
shelf.  Reverberation  in  the  presence  of  the  internal  wave  is  modeled 
numerically.  Calculations  show  that  acoustic  energy  propagating  near  a 
sound  speed  minimum  is  deflected  as  a  high  intensity,  higher  angle  beam 
into  the  bottom,  where  it  is  backscattered  along  the  reciprocal  path.  The 
interaction  of  sound  with  the  internal  wave  is  isolated  in  space,  hence 
resulting  in  a  target-like  clutter,  which  is  found  to  be  greater  than  lOdB 
above  the  mean  reverberation  level. 
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1.  Introduction 

One  major  factor  hindering  the  performance  of  active  systems  is  false  alarms  due  to 
clutter.  Clutter  refers  to  target-like  signals  appearing  in  reverberation.  Clutter  has  been 
attributed  to  strong  backscattering  strength  from,  for  example,  seafloor  features1  and 
fish  schools.2  Another  possible  mechanism  for  clutter,  discussed  in  this  paper,  is  a 
propagation  effect,  with  ordinary  backscattering  strength.  A  sound  speed  feature,  such 
as  a  nonlinear  internal  wave  (NLIW)  with  large  amplitude,  deflects  sound  into  the  bot¬ 
tom  as  a  beam,  causing  a  higher  insonification  at  a  higher  grazing  angle,  both  aspects 
resulting  in  increased  backscatter.  We  demonstrate  that  observed  NLIWs  in  measured 
stratification  can  give  a  greater  than  lOdB  target-like  return  above  the  general  rever¬ 
beration  level.  For  purposes  of  realism,  our  oceanographic  and  acoustic  modeling  is 
based  on  measured  stratification  and  observed  wave  amplitudes  on  the  New  Jersey 
shelf  during  the  SWARM  experiment.3  Ray  tracing  is  performed  that  shows,  qualita¬ 
tively,  the  physical  effect  of  deflection  by  the  internal  wave  of  sound  into  the  bottom 
at  higher  grazing  angle.  A  full  wave  model  then  gives  quantitative  results  for  the 
clutter. 

2.  Swarm  NLIW  model  and  ray  tracing 

Nonlinear  internal  waves  are  very  common  on  the  continental  shelf  where  a  thermo- 
cline  exists.  Several  experiments,  including  SWARM,  combining  oceanography  and 
acoustics  have  been  carried  out  on  the  New  Jersey  shelf.  Furthermore,  this  region  is  of 
continuing  interest  for  additional  experimental  study  in  the  near  future,  therefore,  the 
predictions  presented  here  could  be  investigated  in  future  field  programs. 

The  sound  speed  profile  used  for  modeling  is  derived  from  a  typical  conducti¬ 
vity,  temperature,  and  depth  (CTD)  cast  from  the  SWARM  experiment3  and  is  shown 
in  the  left  panel  of  Fig.  1,  where  a  mixed  layer  extends  down  to  10  m.  A  warm  salty 
layer  near  the  bottom  due  to  the  shelf-break  front  resulted  in  a  minimum  in  the  sound 
speed  at  35  m  depth.  This  is  a  typical  feature  on  the  New  Jersey  coast  and  occurs  in  all 
SWARM  and  Shallow  Water  2006  profiles4  on  the  outer  part  of  the  shelf.  In  this 
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Sound  Speed  (m/s)  Range  (m) 

Fig.  1.  (Color  online)  Sound  speed  profile  (left),  sound  speed  field  with  solitary  wave  at  5  km  in  color  and  super¬ 
imposed  ray  trace  with  grazing  angles  within  ±6  deg.  Almost  all  rays  are  deflected  into  the  bottom  by  the  wave. 


paper,  we  investigate  the  situation  where  a  train  of  nonlinear  internal  waves 
approaches  the  acoustic  source  and  receiver.  In  this  case,  the  leading  wave  usually  has 
the  largest  amplitude,  and  the  sound  speed  profile  is  relatively  range  independent 
before  reaching  that  wave.  The  effect  of  a  leading  wave  of  sufficient  amplitude  is  to 
deflect  the  ducted  sound  near  the  sound  speed  minimum  into  a  beam  that  strikes  the 
bottom,  leading  to  a  target-like  clutter  signal.  The  subsequent  behavior  of  that 
deflected  beam  depends  on  the  details  of  the  following  waves  in  the  train.  Dispersion 
of  the  beam  leads  to  lower  and  more  time  spread  reverberation,  unlikely  to  look  like  a 
target.  Therefore,  we  model  only  the  leading  wave  as  a  solitary  wave,  and  ignore  the 
following  waves.  The  leading  wave  is  almost  a  solitary  wave,  which  we  model  as  a  so¬ 
lution  of  the  fully  nonlinear  Dubriel-Jacotin-Long  (DJL)  equation  (formerly  referred  to 
as  Long’s  equation).5  For  a  given  stratification,  the  DJL  equation  has  a  set  of  one- 
parameter  solutions.  The  amplitude  of  the  wave  used  for  modeling  is  chosen  as  typical 
of  waves  observed  in  SWARM  and  is  shown  as  a  sound  speed  field  in  color  in  Fig.  1 
in  the  online  version.  The  DJL  equation  models  the  two-dimensional  structure  of  a 
NLIW  of  that  amplitude  much  more  accurately  than  does  weakly  nonlinear  equations, 
such  as  the  Korteweg-de Vries  equation.  The  wave  is  assumed  to  be  at  a  5  km  range 
from  the  source  and  receiver  in  order  to  give  quantitative  clutter  to  reverberation. 

To  picture  the  mechanism,  a  set  of  rays  was  launched  from  a  source  at  range 
zero  and  depth  38  m,  close  to  the  sound  speed  minimum.  Figure  1  shows  rays  between 
±6°,  which  stay  in  the  sound  channel  and  contribute  to  the  clutter  mechanism  of  this 
paper,  although  the  detailed  calculations  to  follow  use  an  omnidirectional  source, 
including  propagation  at  the  larger  angles  responsible  for  the  reverberation  at  earlier 
times.  When  the  rays  encounter  the  solitary  wave  they  are  deflected  to  a  higher  angle 
of  about  1 1  deg,  and  leave  the  channel  as  a  beam,  impinging  onto  the  bottom.  This 
beam,  of  higher  grazing  angle  sound,  is  then  backscattered  around  the  reciprocal  path 
to  the  receiver  near  the  source,  arriving  as  a  strong  target-like  clutter  signal.  In  order 
for  this  deflection  of  sound  to  occur,  the  amplitude  of  the  internal  wave  needs  to  be 
large.  Smaller  wave  would  not  change  the  direction  of  the  propagating  sound  out  of 
the  channel.  Ray  tracing  is  unlikely  to  quantitatively  predict  this  phenomenon  for 
frequencies  of  a  few  kilohertz.  Therefore,  we  resort  to  a  full  wave  method  in  Sec.  3  for 
the  quantitative  study;  however,  ray  tracing  still  offers  an  intuitive  picture. 

3.  Reverberation  and  clutter 

For  quantitative  results  on  reverberation,  we  resort  to  a  full  wave  method  to  handle 
the  two-way  propagation,  and  perturbation  theory  for  reverberation.6  The  bottom  is 
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assumed  to  be  a  fluid  half-space  with  sound  speed  of  1650  m/s,  density  of  2000  kg/m3, 
and  an  attenuation  of  0.5  dB/wavelength,  so  seabed  absorption  is  included  in  the  calcu¬ 
lation.  The  acoustic  field  of  an  omnidirectional  point  source  is  expanded  in  unper¬ 
turbed  modes,  and  the  equation  for  their  evolution  including  mode  coupling  by  the 
wave  is  solved  following  Eq.  (29)  of  Ref.  7,  where  we  make  their  forward-scattering 
approximation  by  dropping  the  second  line  and  we  refrain  from  expanding  K2  in  the 
sound  speed  perturbation.  Higher  order  modes  were  included  in  the  calculation  so  the 
high  angle  waves  from  the  point  source  are  properly  included.  Alternatively,  one  could 
use  another  full  wave  approach,  such  as  a  parabolic  equation  method,  to  obtain  the 
same  results.  We  choose  a  point  source  of  3  kHz  at  the  common  choice  of  being  near 
the  sound  axis;  specifically,  our  source  depth  is  38  m.  The  solitary  wave  is  at  a  range 
of  5  km,  and  extends  the  calculation  to  10  km.  Figure  2  shows  the  calculated  horizontal 
flux.  Before  encountering  the  solitary  wave,  most  of  the  sound  energy  propagating  to 
long  distances  is  trapped  in  the  sound  channel  in  the  middle  of  the  water  column.  At 
5  km  range,  the  internal  wave  converts  most  of  the  sound  at  low  modes  (small  angle) 
to  higher  modes  (larger  angle),  in  a  beam  which  interacts  strongly  with  the  seafloor. 
This  result  is  qualitatively  consistent  with  the  ray  analysis  shown  in  Fig.  1.  The 
deflected  beam  is  refracted  by  the  sound  speed  profile,  periodically  hitting  the  bottom. 
At  these  subsequent  bottom  interactions,  the  beam  has  significantly  broadened. 

Target-like  clutter  occurs  when  the  narrow  beam  backscatters  from  the  bot¬ 
tom.  To  calculate  the  reverberation  and  clutter,  the  rough  seafloor  backscatter  model,6 
based  on  perturbation  theory,  is  used.  The  water/sediment  interface  is  assumed  to  have 
small-scale  roughness  with  a  typical  roughness  power  spectrum  as  defined  in  the  Office 
of  Naval  Research  Reverberation  Workshop.8  The  reverberation  vs  time  is  given  in 
Fig.  3,  with  a  100  Hz  bandwidth  centered  on  3  kHz.  At  slightly  before  7  s,  correspond¬ 
ing  to  the  range  where  the  deflected  rays  of  Fig.  1  hit  the  bottom,  a  large  arrival 
appears  more  than  lOdB  above  the  background  reverberation  level  with  a  time  spread 
consistent  with  the  inverse  bandwidth  of  the  transmitted  pulse.  Additional  calculations 
with  twice  the  bandwidth  around  the  same  center  frequency  show  a  2-dB  increase  of 
the  clutter  arrival  while  the  width  remains  to  be  consistent  with  the  inverse  bandwidth. 
The  large  return  is  the  result  of  both  larger  insonification  at  the  place  the  deflected 
beam  strikes  the  bottom  and  the  stronger  backscatter  at  higher  grazing  angles.  Because 
the  deflected  sound  behaves  like  a  narrow  beam  impinging  onto  the  bottom,  the  arrival 
shows  up  like  a  target.  Following  the  large  arrival,  the  reverberation  is  generally  higher 


Range  (m) 

Fig.  2.  (Color  online)  Horizontal  intensity  flux  of  one  way  coupled  mode  solution  at  3000  Hz.  Source  at  38  m 
depth  and  zero  range,  water  column  sound  speed  field  as  in  Fig.  1,  and  sediment  parameters  given  in  text. 
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Time  (s) 

Fig.  3.  Monostatic  reverberation  level  normalized  to  source  level  and  as  a  function  of  time  with  a  point  source 
at  38  m,  5000  m  away  from  the  solitary  wave. 

than  the  level  extrapolated  from  before  the  arrival  of  the  internal  wave.  There  are 
additional  individual  arrivals  due  to  the  repeated  interaction  of  the  deflected  beam 
with  the  bottom.  However,  even  in  the  absence  of  the  rest  of  the  wave  train,  they  are 
more  spread  out  in  time  and  smaller  in  amplitude  than  the  first  peak  because  of  the 
spreading  of  the  deflected  beam. 

4.  Discussion 

Nonlinear  internal  waves  are  a  common  phenomenon  on  the  continental  shelf  when 
strong  stratification  exists.  Under  the  right  conditions,  such  a  wave  induces  a  false  tar¬ 
get  in  an  active  system  due  to  strong  insonification  at  high  angle  of  the  bottom  just 
behind  the  deflecting  wave.  Considering  typical  internal  waves  on  the  shelf  move  at  a 
speed  of  order  1  m/s,  the  induced  false  target  could  be  tracked  as  a  slowly  moving  tar¬ 
get.  When  the  source  is  near  the  channel  minimum,  as  was  shown  in  the  previous 
calculation,  the  effect  is  most  pronounced;  a  source  outside  the  duct  does  not  have  as 
dramatic  an  effect.  To  experimentally  verify  the  effect  presented,  a  monostatic  rever¬ 
beration  system  properly  placed  in  the  duct  aimed  toward  the  oncoming  nonlinear 
internal  waves,  and  a  means  of  monitoring  the  location  of  such  waves,  such  as  radar 
or  visual,  would  suffice. 
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Abstract — During  the  2006  Shallow  Water  (SW06)  experiment, 
simultaneous  measurements  were  made  of  the  sound-speed  field  as 
a  function  of  range  and  depth  associated  with  nonlinear  internal 
waves  and  acoustic  propagation  at  frequencies  of  2-10  kHz  over 
a  1-km  path.  The  internal  waves  were  measured  by  a  towed  con¬ 
ductivity-temperature-depth  (CTD)  chain  to  get  high  resolution. 
These  measurements  were  coordinated  so  that  the  nonlinear  waves 
could  be  interpolated  onto  the  acoustic  path,  allowing  predictions 
of  their  effects  on  the  acoustics.  Using  the  measured  sound-speed 
field,  the  acoustic  arrivals  under  the  influence  of  the  internal  waves 
are  modeled  and  compared  to  data.  The  largest  impact  of  measured 
moderate  amplitude  internal  waves  on  acoustics  is  that  they  alter 
the  arrival  time  of  the  rays  which  turn  at  the  thermocline. 

Index  Terms — Acoustics  propagation,  conductivity-tempera¬ 
ture-depth  (CTD)  chain,  nonlinear  internal  waves. 


I.  Introduction 

THE  2006  Shallow  Water  (SW06)  experiment  was  a 
major  field  effort  conducted  in  July-September  2006  off 
the  coast  of  New  Jersey,  involving  ships,  moorings,  aircraft, 
and  satellite  coverage,  as  well  as  modeling  [1].  Topics  in 
ocean  acoustics,  physical  oceanography,  and  marine  geology 
were  studied.  Oceanographically,  the  shelfbreak  front,  the 
local  eddy  field,  and  the  nonlinear  internal  wave  field  were 
of  interest.  Acoustically,  SW06  included  both  low-frequency 
(50-1600  Hz)  and  midfrequency  (1600-20000  Hz)  compo¬ 
nents,  and  examined  issues  in  forward  propagation,  scattering, 
and  inverse  theory.  For  acoustics  transmissions,  both  moored 
and  shipboard  sources  and  receivers  were  used.  The  low-fre¬ 
quency  component  concentrated  on  two  issues:  geoacoustic 
inversion  and  azimuthal  dependence  of  shallow-water  sound 
propagation,  specifically  its  strong  dependence  upon  the  angle 
between  the  acoustic  path  and  the  propagation  direction  of 
the  nonlinear  internal  wave  field.  The  midfrequency  efforts 
of  SW06  were  concentrated  on  receiving  arrays  moored  in 
the  central  area  where  the  water  depth  is  about  80  m.  In  this 
central  area,  extensive  environmental  measurements  were  made 
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to  support  midfrequency  acoustics.  Measurements  include 
short-range  acoustics  interaction  with  the  bottom  and  surface, 
short-range  propagation  through  internal  waves,  and  longer 
range  propagation  to  10  km  along  the  same  track  with  multiple 
source  depths,  frequency  bands,  and  times.  The  location  of 
the  nominal  center  of  the  experiment  was  39°  N,  73°  W.  In 
this  paper,  we  choose  this  point  to  be  (x.  y)  =  (0, 0)  with  x 
pointing  to  the  east  and  y  pointing  to  the  north.  One  of  the  major 
themes  of  the  program  was  highly  nonlinear  internal  waves 
(NLIWs).  While  oceanographic  studies  were  concerned  with 
the  creation,  propagation,  and  decay  of  those  waves,  acoustic 
studies  were  concerned  with  their  effects  on  propagation.  This 
paper  reports  a  study  of  sound  propagation  through  NLIWs, 
which  were  concurrently  measured  to  enable  deterministic 
model/data  comparisons.  This  study  was  originally  planned  for 
August  1 1 ,  but  it  turned  out  (very  unusually)  that  no  mode  1 
NLIWs  were  observed  near  the  experimental  site  that  day.  It 
was  rescheduled,  and  performed  on  August  13. 

The  impact  of  internal  waves  on  sound  propagation  in  coastal 
areas  has  been  studied  quite  extensively  in  the  past  two  decades. 
For  examples,  SWARM  [2]  on  the  New  Jersey  shelf  north  of  the 
SW06  site,  the  Slope  to  Shelf  Primer  [3],  and  the  2001  Asian 
Sea  International  Acoustics  Experiment  (ASIAEX)  in  the  South 
China  Sea  [4] .  These  programs  have  concentrated  on  propaga¬ 
tion  of  low-frequency  sound  (<1  kHz)  [5].  Mode  conversion 
and  acoustic  propagation  in  directions  other  than  the  NLIW 
propagation  direction  dominate  these  low-frequency  studies. 

One  part  of  SW06  was  designed  to  provide  sufficient  oceano¬ 
graphic  measurements  of  NLIWs  so  that  the  2-D  sound-speed 
fields  (range  and  depth)  in  the  acoustic  path  at  the  time  of 
sound  transmissions  could  accommodate  deterministic  mod¬ 
eling  of  the  effects  of  NLIWs  on  sound  propagation.  The 
instrument  used  for  measuring  these  NLIWs  was  a  towed 
conductivity-temperature-depth  (CTD)  chain,  which  provides 
a  2-D  sound-speed  field  over  time.  The  towed  CTD  chain  data 
are  augmented  by  conventional  CTD  casts. 

Under  the  usual  summer  conditions  with  a  thinner  warm 
upper  mixed  layer,  a  strong  thermocline,  and  a  thicker,  less 
stratified  lower  layer,  NLIWs  are  waves  of  depression.  In 
SW06,  the  waves  had  the  common  property  of  coming  in  wave 
trains  with  larger  waves  tending  to  occur  near  the  front  of  the 
train  [1,  Fig.  2].  NLIWs  in  the  experiment  moved  with  a  speed 
between  0.5  and  1.0  m/s,  and  had  horizontal  widths  perpen¬ 
dicular  to  the  propagation  direction  of  the  NLIW  on  the  order 
of  100  m.  These  result  in  variations  of  sound  propagation  over 
time.  To  model  such  variations,  the  sound- speed  field  along  the 
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acoustic  path  must  be  known  as  a  function  of  time.  The  speed 
and  spatial  structure  of  the  NLIWs  make  it  impractical  to  move 
a  single  CTD  vertically  through  the  water  fast  enough  to  resolve 
a  wave’s  structure  well.  Multiple  CTDs  positioned  vertically 
are  needed  for  good  resolution.  They  can  either  be  mounted  on 
a  mooring  or  lowered  from  a  ship  or  float.  Viewed  from  the 
reference  frame  moving  with  the  wave,  the  wave’s  evolution 
over  time  is  slow.  Hence,  a  towed  chain  has  the  advantage  over 
a  mooring  because  the  slow  wave  evolution  can  be  measured 
(and  was  in  SW06),  in  addition  to  the  rapid  time  dependence 
with  the  wave  speed. 

The  planning  of  the  acoustics  experiment  was  based  on  his¬ 
torical  NLIW  data  on  the  New  Jersey  shelf  in  the  summer.  The 
summer  sound- speed  profiles  of  the  area  typically  have  a  rel¬ 
atively  thin  mixed  layer  near  the  surface  with  a  high  sound 
speed,  a  thin  thermocline  to  reach  a  minimum  sound  speed,  and 
a  thicker  region  of  less  variable  sound  speed  extending  to  the 
bottom.  At  the  time  and  place  of  concern  in  this  paper,  as  fre¬ 
quently  happens  in  this  area,  a  warm,  salty  bottom  layer  intruded 
from  off  the  shelf,  making  the  sound-speed  gradient  near  the 
bottom  somewhat  higher  than  with  an  adiabatic  profile.  The  re¬ 
sult  is  a  sound  channel  with  higher  speed  both  near  the  surface 
and  near  the  bottom.  Such  a  channel  favors  propagation  and  re¬ 
duces  the  influence  of  the  surface  and  bottom,  especially  if  the 
sound  source  is  placed  near  the  sound  axis.  At  the  experimental 
site,  the  water  depth  is  80  m.  The  sound  axes  varied  in  depth 
from  30  to  40  m  during  the  experiment. 

Acoustic  intensity  fluctuations  are  largely  due  to  spatially 
small  sound-speed  perturbations,  and  a  deterministic  interpo¬ 
lation  of  such  perturbations  is  unlikely  to  give  correct  results. 
Rather,  we  concentrate  on  phenomena  sensitive  to  large  spatial 
scales,  which  are  associated  with  the  travel  time  of  acoustic  ar¬ 
rivals. 

Historical  oceanographic  data  suggested  that  NLIWs  in  the 
area  come  in  trains  [2]  with  spacing  on  the  order  of  100  m  and 
they  propagate  from  the  southeast  to  northwest  generally  in  the 
direction  along  a  300°  compass  bearing. 

The  sound  source  and  receiver  were  separated  by  1000  m 
along  the  300°  compass  direction.  The  CTD  chain,  towed  by  the 
R/V Endeavor,  circled  close  to  the  acoustic  path,  typically  200  m 
from  the  path,  while  midfrequency  sound  propagation  data  were 
collected.  The  tow  track  is  shown  in  Fig.  1  together  with  the 
source  and  receiving  array  positions.  The  source  was  suspended 
from  the  stern  of  the  R/V  Knorr  that  was  dynamically  posi¬ 
tioned  at  39.0203°  N,  73.0277°  W,  facing  due  north,  putting  the 
source  at  (#,  y)  =  (—2.3999,2.2133)  km,  40  m  due  south  of  the 
global  positioning  system  (GPS)  position.  The  receiving  array, 
moored  at  39.0245°  N,  73.0377°  W,  is  about  1000  m  from  the 
source,  at  a  bearing  of  300°.  During  this  joint  acoustic/oceano¬ 
graphic  measurement  period,  NLIWs  crossed  the  acoustic  path, 
allowing  the  study  of  their  impact  on  sound  propagation. 

In  the  following,  we  first  present  the  oceanic  and  acoustic 
data  in  Section  II,  and  then  in  Section  III,  we  model  the  oceanic 
data  to  interpolate  the  sound-speed  field  onto  the  acoustic  path. 
In  Section  IV,  we  model  the  acoustic  propagation  through  the 
sound-speed  field.  In  Section  V,  we  attempt  to  interpret  the  sys¬ 
tematic  features  of  acoustic  arrivals  using  ray  trace  ideas.  We 
finish  with  a  discussion  in  Section  VI. 


E-W  (km) 

Fig.  1.  Positions  of  the  acoustic  and  oceanographic  measurements.  Distances 
are  relative  to  longitude  73°W  and  from  39°N.  The  R/V  Endeavor  track, 
16:42:00Z  to  17:15:00Z  August  13,  circling  around  the  acoustic  source  and 
receiving  array.  Thin  solid  line:  leg  1,  16:42:00Z  to  16:53: 00Z  ;  thick  solid  line: 
leg  2,  16:53:00Z  to  17:04:00Z;  and  dashed  line:  leg  3,  17:04:00Z  to  17:15:00Z. 
Encounters  with  the  center  of  the  first  wave  are  indicated  by  the  diamond, 
circle,  and  triangle  for  legs  1,  2,  and  3,  respectively. 

II.  Experimental  Data 

A.  Oceanic  Data 

During  the  acoustic  transmissions,  the  density  and  sound 
speed  were  measured  by  a  CTD  chain  towed  by  the  R/V  En¬ 
deavor. 

The  towed  CTD  chain  system  comprises  a  set  of  individual 
CTD  fins  arranged  on  a  single  cable  that  is  towed  behind  a  ship. 
The  chain  is  deployed  from  a  2-m  diameter  divided  drum.  To 
facilitate  recovery,  the  drum  is  optimally  situated  so  that  the 
fins  come  off  the  bottom  of  the  drum,  allowing  the  fins  to  orient 
themselves  properly  with  the  drum  during  recovery.  The  chain 
threads  through  a  specially  designed  block  suspended  from  the 
A-frame.  To  maximize  the  chain  depth  while  towing,  a  large 
V-wing  depressor  is  attached  to  the  end  of  the  chain  cable. 

The  CTD  fins  are  molded  rubber  approximately  190  mm  long 
by  100  mm  tall  with  a  rounded  body  that  tapers  from  50  to 
30  mm.  The  insulated  cable  acts  as  both  strength  member  for 
towing  and  as  power/communication  line  for  data  collection.  In¬ 
dividual  sensors  are  threaded  onto  the  tow  cable  through  a  hole 
parallel  to  the  100-mm  side,  at  the  wide  end  of  the  body.  The 
sensors  are  free  to  rotate  around  the  cable  so  that  they  self-align 
to  decrease  drag  while  being  towed.  The  sensors  are  distributed 
along  the  wire  according  to  the  number  of  sensors  available  and 
the  desired  sampling  in  depth.  They  are  held  in  place  using  split 
collars  clamped  to  the  cable  5  mm  above  and  below  each  sensor. 
Each  CTD  fin  contains  sensors  to  measure  pressure  (P),  temper¬ 
ature  (T),  and  conductivity  (C).  Power  and  data  transmission  are 
achieved  through  inductive  coupling  of  each  sensor  with  the  tow 
cable. 

As  configured  on  August  13,  48  sensors  were  spaced  about 
1  m  vertically  (1.5  m  along  the  slant  of  the  chain),  and  the  chain 
was  towed  at  about  6  kn.  GPS  positions  and  time  and  the  mea¬ 
sured  pressure  converted  to  depth  are  used  for  the  determination 
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Fig.  2.  Fields  obtained  from  the  towed  CTD  chain,  (a)  Density  — 1000  kg/m3  on  the  first  shoreward  leg  (leg  1);  (b)  sound  speed  corresponding  to  (a);  (c)  sound 
speed  from  seaward  leg  (leg  2);  (d)  sound  speed  from  the  second  shoreward  leg  (leg  3).  The  source  and  receiver  locations  relative  to  the  waves  are  also  shown  for 
each  leg.  The  density  and  sound  speed  show  the  same  features,  except  that  the  density  shows  calibration  inaccuracies.  Similar  relationships,  not  shown,  hold  on 
the  other  two  legs.  We  use  sound  speed  as  a  surrogate  for  density.  The  same  two  nonlinear  internal  waves  show  up  on  all  three  legs. 


of  the  actual  sensor  positioning.  The  measurements  dedicated 
to  acoustics  lasted  for  2.5  hours,  during  which  the  CTD  chain 
was  towed  in  seven  circuits  around  the  propagation  path.  At  this 
time,  there  were  two  simultaneous  trains  of  NLIWs  along  the 
acoustic  path.  These  two  trains  merged  into  one  train  along  a 
line  that  passed  a  short  distance  southwest  of  the  circuits,  and  the 
merge  was  clearly  visible  from  the  RJV  Endeavor  from  which 
the  CTD  chain  was  deployed.  After  the  acoustics  support  phase, 
the  RJV  Endeavor  proceeded  northwest  to  the  merge  region  of 
the  fronts  of  the  trains,  to  study  the  merging  process. 

One  and  a  half  circuits  (which  we  call  three  legs  shown  in 
Fig.  1)  in  the  middle  of  this  time  period  are  selected  for  analysis, 
as  the  large  waves  at  the  front  of  the  packets  passed  the  acoustic 
path  during  this  time. 

The  legs  extended  some  distance  beyond  the  source  and  the 
receiver.  This  was  done  so  that  the  chain  would  have  a  con¬ 
stant  shape  and  be  directly  behind  the  RJV  Endeavor  when  the 
track  was  along  the  acoustic  path.  The  measured  density  and 
sound  speed  on  the  first  leg  are  shown  in  Fig.  2(a)  and  (b).  The 
sound  speed  is  nearly  a  unique  function  of  the  density;  very 
little  “spice”  variability  occurs.  This  is  true  of  the  other  two 
legs  as  well,  and  only  the  measured  sound  speed  is  shown  in 
Fig.  2(c)  and  (d).  By  comparing  Fig.  2(c)  and  (d)  to  their  corre¬ 
sponding  densities,  we  find  that  they  also  have  very  little  “spice” 
variability  and  the  sound-speed-density  relation  is  the  same  as 
that  of  the  first  leg.  Two  large  waves  are  seen  in  all  three  legs. 
The  spacing  between  these  waves  is  the  same  in  the  first  and 


third  legs,  but  is  less  in  the  second  leg,  even  after  accounting  for 
the  relative  motion  of  the  ship  and  the  waves  on  all  legs.  The 
shape  of  the  second  wave  is  also  different  in  the  second  leg.  The 
second  leg  is  400  m  to  the  southwest  of  legs  1  and  3,  closer  to 
the  merging  region  of  the  wave  trains,  so  these  differences  are 
along-crest  variability  of  the  waves. 

One  profile  taken  at  15:30:00Z  on  August  13  with  the  ship’s 
CTD  on  the  RJV Knorr  is  important  in  our  analysis.  This  profile 
was  taken  at  the  approximate  location  of  the  acoustic  source  in 
our  experiment,  and  is  shown  in  Fig.  3(a).  For  reasons  of  cost, 
CTD  units  on  the  towed  chain  are  not  as  high  quality  as  a  ship’s 
CTD.  Accuracy,  reliability,  and  calibrations  (especially  relative 
calibrations  between  different  fins)  are  issues  that  have  not  been 
adequately  resolved  for  detailed  quantitative  use.  Therefore,  the 
ship’s  CTD  profile  was  used  to  complement  the  towed  chain 
data  in  obtaining  the  sound-speed  field  for  acoustics  simula¬ 
tions.  The  method  of  combining  the  data  is  described  below. 


B.  Acoustic  Data 

The  acoustic  range  is  nominally  1000  m,  the  source  was  po¬ 
sitioned  at  30  m  depth,  and  the  receiver  used  in  this  analysis 
was  25  m  deep.  The  total  water  depth  is  80  m  at  the  experi¬ 
mental  site.  The  acoustic  source  was  deployed  from  the  stern 
of  the  RJV  Knorr ,  which  was  dynamically  positioned  at  a  fixed 
point.  The  acoustic  path  was  at  a  bearing  of  300° .  Because  of  the 
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Fig.  3.  (a)  Sound  speed  from  a  CTD  profile  taken  at  15:30:00Z  on  August  13  near  the  source  position  as  given  in  Fig.  1.  The  arrow  points  to  1510  m/s,  which 
is  used  to  define  the  thermocline  depth.  Small  errors  in  this  sound-speed  value  in  the  towed  chain  measurement  would  make  little  difference  to  the  thermocline 
depth,  (b)  Rays  between  the  source  and  the  receiver  using  a  range-independent  sound  speed  extracted  from  this  CTD  profile.  The  rays  are  named  bottom-turn  (BT), 
turn-bottom  (TB),  turn-bottom-turn  (TBT),  and  surface-bottom-surface  (SBS). 


source  deployment  method,  there  is  an  uncertainty  of  0(1  m)  in 
the  acoustic  range. 

The  transmitted  acoustic  signal  spanned  the  frequencies 
1.5  and  10.5  kHz.  A  set  of  direct-path  propagation  data 
was  collected  at  the  50-m  range  as  calibration,  so  the  trans¬ 
mitted  waveforms  were  known.  In  the  frequency  band,  the 
source  has  a  10-dB  notch  near  6  kHz.  However,  because  the 
signal-to-noise  ratio  (SNR)  was  high  at  the  1-km  range,  data 
from  the  full  band  can  be  utilized.  The  pulse  was  compressed 
by  dividing  the  Fourier  transform  of  the  received  pulse  by 
the  Fourier  transform  of  the  sent  pulse  and  multiplying  by 
sin(7 r((/  —  1.5  kHz)/9  kHz))  between  1.5  and  10.5  kHz.  The 
time  resolution  is  estimated  as  the  full  width  at  half  maximum 
intensity  of  this  compressed  analytic  signal,  which  is  0.13  ms. 
The  sidelobes  are  down  by  23  dB. 

Acoustic  transmission  started  at  14: 14:00Z  and  ended 
18:49:00Z  on  August  13  with  a  nominal  14-s  ping  repeti¬ 
tion  rate,  resulting  in  1100  pings  recorded.  The  uncertainty 
of  the  source  position  excluded  the  possibility  of  coherent 
ping-to-ping  processing  and  motivated  use  of  the  highest  arrival 
of  each  ping  as  reference  time  zero  to  line  up  all  the  pings. 
The  highest  arrival  almost  always  occurred  during  the  arrival 
of  sound  propagating  near  the  minimum  sound-speed  depth. 
Internal  waves  have  little  effect  on  the  sound  speed  near  its 
minimum,  making  this  a  good  reference.  We  define  the  reduced 
travel  time  as  the  actual  travel  time  shifted  by  the  amount 
corresponding  to  the  reference  time.  In  Fig.  4,  acoustic  arrival 
intensity  for  the  25 -m  receiver  is  plotted  against  geotime  on  the 
vertical  axis  and  the  reduced  travel  time  on  the  horizontal  axis 
for  all  the  transmissions. 
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Fig.  4.  Acoustic  intensity  level  during  this  experiment.  The  level  is  relative  to 
the  highest  intensity  observed  during  this  period.  Six  arrivals  are  labeled  near 
the  top  of  the  figure.  Arrival  #1  is  a  faint  arrival  that  has  propagated  through  the 
surface  mixed  layer,  arrival  #2  is  the  direct  arrival  used  to  define  the  reduced 
travel  time,  and  arrivals  #3-#6  are  the  arrivals  associated  (by  travel  time)  with 
the  rays  of  Fig.  3(b).  The  correspondence  is  #3  =  BT,  #4  =  TB,  #5  =  TBT,  and 
#6  =  SBS.  The  oscillations  in  arrival  time  around  16:30:00Z,  most  clearly  seen 
in  arrival  #5,  are  the  primary  concern  of  our  modeling. 


In  the  30-ms  time  window  shown  in  Fig.  4,  six  separated 
arrivals  are  identified  and  labeled.  Between  geotime  14:20:00Z 
and  16:00:00Z,  before  the  waves  are  present,  the  arrival  struc¬ 
ture  is  relatively  stable:  arrival  #1  being  very  weak  at  about 
—  15  ms,  followed  by  arrival  #2  with  the  highest  amplitude 
defining  0  ms,  in  turn  followed  by  arrivals  #  3,  #4,  and  #5, 
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between  2  and  5.5  ms,  ending  by  a  weaker  arrival  #6  at  about 
7.0  ms. 

Ray  tracing  was  performed  using  a  sound-speed  profile  from 
a  CTD  cast  from  the  R/V  Knorr  at  15:30:00Z  geotime  as  given 
in  Fig.  3(a)  to  identify  the  different  arrivals.  The  rays  associ¬ 
ated  with  each  of  the  labeled  arrivals  in  Fig.  4  based  on  travel 
time  are  as  follows.  Arrival  #1,  a  faint  arrival,  is  a  ray  leaving 
the  source  upward,  traveling  near  the  surface,  and  reaching  the 
receiver  going  downward.  Arrival  #2,  with  the  highest  ampli¬ 
tude,  is  a  superposition  of  two  kind  of  rays;  the  first  kind  is  a 
group  of  axial  arrivals  traveling  close  to  the  sound  axis  at  30-m 
depth  without  reaching  the  high  gradient  part  of  the  thermo- 
cline  and  the  second  kind  is  a  ray  labeled  B  in  the  convention  of 
Fig.  3,  leaving  the  source  downward,  bouncing  off  the  bottom 
once  then  reaching  the  receiver  upward.  The  two  kinds  of  ar¬ 
rivals  reach  the  receiver  spread  out  in  time  over  1  ms,  hence 
are  not  resolved.  The  next  four  labeled  arrivals  are  identified 
with  the  rays  shown  in  Fig.  3(b).  Arrival  #3  starts  downward, 
bounces  once  off  the  bottom,  then  turns  downward  in  the  ther- 
mocline,  reaching  the  receiver  going  downward  (BT).  Arrival 
#4  can  be  considered  a  reciprocal  of  arrival  #3,  which  starts  up¬ 
ward,  turns  downward  in  the  thermocline,  then  bounces  at  the 
bottom,  reaching  the  receiver  going  upward  (TB).  Arrival  #5  is 
a  ray  that  has  two  upper  turning  points  and  one  bottom  bounce 
between  them  (TBT).  Arrival  #6  is  a  ray  bouncing  twice  at  the 
sea  surface  and  once  at  the  bottom  (SBS). 

The  arrivals  in  Fig.  4  have  a  changing  structure  versus  geo¬ 
time.  They  tend  to  get  closer  to  one  another  starting  from  about 
16:30:00Z.  Notably,  three  of  the  arrival  times  have  oscillatory 
patterns.  The  “TBT”  arrival  has  four  fluctuations  to  earlier  ar¬ 
rival  times,  whereas  the  “BT”  and  “TB”  arrivals  have  two  such 
fluctuations  in  each.  However,  the  arrival  “SBS,”  though  it  has  a 
similar  ray  path  to  “TBT,”  does  not  show  apparent  fluctuations 
in  its  arrival  time.  In  addition,  the  overall  spread  in  arrival  times 
is  reduced  by  about  1  ms  after  16:30:00Z,  with  the  reduction 
occurring  between  arrivals  #4  and  #5. 


III.  Environmental  Modeling 

The  first  step  of  modeling  the  sound-speed  field  is  to  interpo¬ 
late  the  CTD  chain  data  in  space  and  time  to  estimate  the  posi¬ 
tions  of  the  waves  along  the  acoustic  path.  The  wave  speeds, 
their  directions,  and  their  time  crossing  a  reference  position 
[which  we  extrapolate  to  (x,  y)  =  (0,0)]  need  to  be  determined. 
The  observed  spacing  between  the  waves  on  the  two  shoreward 
legs  is  nearly  the  same  (and  the  ship  speed  was  the  same),  so  the 
two  waves  are  traveling  at  about  the  same  speed.  The  wave  speed 
in  the  direction  of  the  acoustic  path  is  primarily  determined  by 
the  position  and  time  spacing  between  the  observations  of  the 
same  wave  on  consecutive  shoreward  passes.  The  wave  widths, 
and  the  distance  between  them,  differ  from  the  raw  observation 
because  the  ship  speed  is  only  several  times  the  wave  speed.  The 
wave  propagation  directions  also  depend  on  the  measurements 
on  the  seaward  leg.  Rather  than  a  sequential  determination  of  the 
parameters,  we  fit  all  simultaneously  from  the  data  on  all  three 
legs  and  interpolate  to  the  acoustic  path,  avoiding  assumptions 
such  as  that  the  tow  path  be  exactly  parallel  to  the  acoustic  path. 
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There  is  a  complex  structure  at  smaller  scales  than  the  two 
waves,  and  it  is  unlikely  that  the  smaller  structure  is  traveling  in 
the  same  direction  with  the  same  speed  as  the  large  waves.  Thus, 
we  only  want  to  interpolate  the  large  scale  features.  We  define 
the  center  of  the  thermocline  by  sound  speed  =1510  m/s,  as 
indicated  in  Fig.  3(a).  We  refer  to  the  depth  at  which  this  is  the 
sound  speed  as  the  “thermocline  depth.”  For  every  towed  chain 
measurement  (every  2  s),  we  estimate  the  thermocline  depth.  We 
then  lowpass  filter  the  resulting  time  series  of  depths  to  elimi¬ 
nate  the  small  waves.  This  lowpass  filter  also  eliminates  about 
0.30  m  from  the  deepest  excursions  of  the  large  wave  (where 
the  smallest  horizontal  scales  occur),  and  this  extra  depth  was 
replaced  at  the  end  of  the  processing. 

In  each  case,  the  lowpass  curve  intersects  the  raw  curve  near 
the  places  where  the  curves  are  2  m  above  the  deepest  point  on 
the  lowpassed  first  wave.  For  each  leg,  these  points  are  taken 
to  define  the  front  and  the  back  of  the  two  waves,  except  that 
the  back  of  the  second  wave  can  only  be  determined  on  leg  3. 
On  leg  1,  the  back  was  outside  the  tow  path,  while  on  leg  2,  the 
thermocline  stayed  deep  for  some  distance  behind  the  second 
wave.  It  is  assumed  that  the  front  and  the  back  of  the  first  wave 
and  the  front  of  the  second  wave  move  at  constant  speed  and 
they  form  straight  lines.  The  three  legs  provide  exactly  enough 
information  to  solve  for  each  speed,  direction,  and  position  of 
each  wave  feature  examined.  Table  I  summarizes  the  solutions. 
The  propagation  directions  are  a  little  more  northward  than  the 
300°  assumed  in  pre-experiment  modeling,  as  were  most  of  the 
waves  observed  in  SW06.  The  intersections  of  the  lines  with  the 
acoustic  path  are  the  modeled  positions  of  those  features  along 
that  path.  The  speeds  as  observed  along  the  path  of  all  three 
points  are  nearly  the  same,  so  we  use  their  average,  0.6447  m/s, 
fixing  the  positions  at  the  time  16:59:00Z,  when  the  waves  were 
near  the  center  of  the  acoustic  path.  The  maximum  difference 
of  any  feature  on  the  acoustic  path  between  using  the  individual 
speeds  and  using  the  average  speed  was  only  15  m.  Finally,  the 
assumed  back  of  the  second  wave  was  chosen  so  that  the  ratio  of 
widths  of  the  two  waves  was  the  same  as  the  ratio  observed  on 
the  third  leg,  very  nearly  parallel  to  the  acoustic  path.  The  center 
of  each  wave  is  placed  half  way  between  its  front  and  back. 

In  addition  to  the  two  waves  shown  in  Fig.  2,  the  towed  chain 
data  show  a  change  in  the  thermocline  depth,  which  changed 
from  about  15  m  before  the  waves  to  about  17  m  after  the  wave. 
These  estimates  ignore  smaller  waves  that  appear  in  the  data,  but 
cannot  be  connected  between  the  different  legs.  This  change  of 
thermocline  depth  is  the  solibore  aspect  of  a  train  of  NLIWs,  as 
discussed  by  Henyey  and  Hoering  [6].  Between  the  two  waves, 
the  lowpass  curve  shows  the  thermocline  depth  at  about  16  m, 
half  way  through  the  bore.  The  thermocline  depth  at  the  centers 
of  the  waves  is  linearly  interpolated  between  the  lowpass  values 
on  the  three  legs,  with  an  extra  depth  of  0.30  m  added.  This  extra 
thermocline  depth  makes  its  integral  between  the  front  and  the 
back  of  the  first  wave  very  nearly  equal  to  that  of  the  raw  curve 
on  each  of  the  three  legs.  Thus,  we  have  seven  points  specified 
on  the  thermocline  depth,  as  well  as  the  two  asymptotes.  With 
the  seven  points  on  the  curve  identified,  plus  the  two  asymp¬ 
totes,  the  remainder  of  the  curve  is  obtained  by  Hermite  inter¬ 
polation  using  the  algorithm  of  Fritsch  and  Carlson  [7] .  This 
method  places  the  local  extrema  at  the  centers  of  the  waves  and 
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TABLE  I 

Parameters  of  Measured  NLIW  Features  Inferred  From  Their  Observation  by  the  Towed  CTD  Chain  on  Three  Legs.  The  Last  Row  Is  Speed  as 
Observed  Along  the  Acoustic  Path  at  a  Bearing  of  300° .  They  Are  Deduced  From  the  Speed  and  Wave  Bearing  in  Rows  1  and  2.  The  Width  of 
the  Second  Wave  Between  Front  and  Back  Is  Assumed  to  Have  the  Same  Ratio  to  the  Width  of  the  First  Wave  as  Observed  on  Leg  3 


Front  of  first 

wave 

Back  of  first 

wave 

Front  of  second 

wave 

Speed 

0.587  m/s 

0.606  m/s 

0.564  m/s 

Bearing 

324.3° 

316.8° 

331.4° 

Extrapolated  time 
at  73°W,  39°N 

15:14:26Z 

15:19:57Z 

15:23:01Z 

Speed  along 
acoustic  path 

0.643  m/s 

0.632  m/s 

0.658  m/s 

Fig.  5.  Model  of  the  thermocline  depth  (at  16:59:00Z)  defined  by  sound  speed 
=  1510m/sasa  function  of  range  obtained  by  interpolating  the  towed  chain  data 
onto  the  acoustic  path.  Circles:  interpolated  points  using  the  tow  chain  data; 
crosses:  auxiliary  points  selected  to  make  the  curve  smooth.  At  the  reference 
time  16:59:00Z,  the  source  is  at  range  0  m,  and  the  receiver  is  at  1000  m. 


Fig.  6.  True  displacement  (solid)  at  the  center  of  a  large  NLIW  compared  to  the 
linear  displacement  (dashed)  (a)  in  the  Eulerian  reference  frame  and  (b)  in  the 
Lagrangian  reference  frame.  The  Lagrangian  reference  frame  is  clearly  a  much 
better  choice  in  which  to  use  the  linear  mode  displacement. 


at  the  16-m  point  between  the  waves,  as  we  intend.  To  ensure  a 
smooth  transition  from  the  waves  to  the  asymptotes,  five  auxil¬ 
iary  points  were  added  in.  Fig.  5  shows  the  result  of  this  inter¬ 
polation. 

Given  the  thermocline  depth  as  a  function  of  the  horizontal 
position  and  time,  we  complete  the  specification  of  the  sound 
speed  using  the  data  from  the  CTD  profile  at  15:30:00Z.  We 
assume  that  the  vertical  displacements  from  that  profile  are 
linear  mode  1  in  a  vertically  Lagrangian  coordinate  system, 
i.e.,  c(x,  zo  +  a(x)<pi(zo))  =  co(zo),  where  is  the  mode 
1  displacement  eigenfunction  and  Co(zo)  is  the  sound-speed 
profile,  both  from  the  CTD  profile  at  15:30:00Z.  This  assump¬ 
tion  is  rather  accurate  for  NLIWs  when  the  profile  has  a  sharp 
thermocline  with  relatively  little  stratification  above  and  below 
[13].  Such  a  quasi-2  layer  profile  describes  well  the  conditions 
during  our  experiment.  Following  [13],  Fig.  6  compares  the 
linear  mode  1  eigenfunction  with  the  displacement  at  the  center 


of  a  large  NLIW  in  two  coordinate  systems,  demonstrating  the 
superiority  of  using  the  linear  displacement  in  the  vertically 
Lagrangian  coordinate  system  over  that  in  the  Eulerian  coordi¬ 
nate  system.  A  vertically  Lagrangian  coordinate  system  is  also 
appropriate  in  other  internal  wave  cases  [8].  The  strength  of  the 
linear  wave  function  at  each  value  of  x  —  vt  is  chosen  to  place 
the  thermocline  depth  on  our  interpolated  curve.  This  model 
sound- speed  field  at  16:59:00Z  is  shown  in  Fig.  7. 

IV.  Acoustic  Modeling 

The  acoustic  data  show  that  when  an  NLIW  train  goes 
through  the  acoustic  path,  different  arrivals  experience  different 
degrees  of  variations.  The  most  apparent  variation  is  that 
some  of  the  arrivals  show  oscillatory  arrival  time  on  the  order 
of  1  ms,  whereas  others  do  not.  To  understand  the  observed 
change,  both  broadband  Collins-type  “parabolic  equation”  (PE) 
[9]  and  ray  tracing  techniques  are  used  to  simulate  acoustic 
propagations  in  this  changing  environment.  While  ray  tracing 
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Fig.  7.  Sound-speed  field  generated  using  the  thermocline  depth  shown  in 
Fig.  5.  The  vertical  dependence  at  each  point  is  assumed  to  be  given  by  a 
vertically  Lagrangian  mode  1  displacement  of  the  profile  shown  in  Fig.  3(a). 
At  the  reference  time  16:59:00Z,  the  source  is  at  range  0  m,  and  the  receiver 
is  at  1000  m. 


Fig.  8.  Comparison  of  acoustic  data  with  a  range-independent  broadband  PE 
simulation  based  on  the  profile  in  Fig.  3.  The  arrivals  are  well  fit  at  the  time  of 
the  profile;  no  NLIWs  were  present  at  that  time.  Even  the  faint  arrival  at  — 15  ms 
is  fit. 


has  the  advantage  of  providing  intuitive  physical  interpretations 
of  observations,  the  PE  method  applies  to  situations  in  which 
diffraction  is  important,  where  ray  tracing  might  be  inaccurate 
or  inapplicable.  A  disadvantage  of  the  PE  simulation  is  that  it 
does  not  easily  associate  features  of  the  sound-speed  field  with 
features  of  the  acoustic  arrivals,  as  ray  tracing  does. 

A.  PE  Modeling  Results 

The  purpose  of  the  PE  simulation  is  to  verify  if  the  overall 
acoustic  arrivals  identified  by  ray  tracing  [Fig.  3(b)]  can  be  re¬ 
produced  both  in  their  arrival  structure  and  in  their  relative  in¬ 
tensity.  The  ocean  environmental  input  to  the  PE  simulations 
is  a  range-independent  model  using  the  CTD  profile  [Fig.  3(a)] 
sound  speed.  The  ocean  surface  is  assumed  flat  and  a  flat  fluid 
half-space  bottom  was  assumed  where  the  geoacoustic  parame¬ 
ters  are  as  follows:  sound  speed  1620  m/s,  density  1850  kg/m3, 
and  absorption  0.5  dB/A  [10]. 

Broadband  pulses  were  generated  by  Fourier  synthesizing 
monochromatic  PE  runs.  In  the  simulation,  the  starting  fre¬ 
quency  is  1.5  kHz,  the  ending  frequency  is  10.5  kHz,  and  the 
frequency  increment  is  2.5  Hz.  These  parameters  result  in  a 
400-ms  time  window,  long  enough  to  encompass  all  desired 
arrivals  without  wrap  around.  In  Fig.  8,  the  envelope  of  the 
synthesized  pulse  is  shown  as  the  white  line  and  compared  with 
acoustic  data  recorded  around  the  time  the  CTD  profile  was 
taken.  The  simulated  waveform  matches  the  acoustic  data  in 
terms  of  number  of  arrivals,  arrival  times,  and  relative  strengths, 
including  the  very  weak  early  arrival  near  —15  ms.  The  PE 
simulation  confirms  that  the  acoustic  arrivals  are  well  modeled 
when  in  situ  oceanographic  and  geophysical  input  is  used.  To 
capture  and  interpret  the  variation  of  the  acoustic  arrivals  when 
NLIWs  are  present,  we  rely  on  ray  tracing  simulations. 

The  broadband  PE  method  shows  good  agreement  with  the 
acoustic  data  at  the  time  the  CTD  profile  was  taken.  Presum¬ 
ably,  it  would  work  well  at  other  times  to  the  extent  that  the 
sound-speed  field  was  properly  modeled,  including  the  fluctua¬ 
tions  in  the  arrival  times.  However,  there  is  no  easy  way  from 
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Fig.  9.  Comparison  of  acoustic  data  with  ray  tracing  results  for  arrivals  (from 
left  to  right  as  indicated  by  the  white  curves):  TB,  BT,  TBT,  and  SBS  between 
16:24:00Z  and  17:26:00Z.  Only  one  overall  reference  time  was  adjusted,  as  the 
axial  arrival  was  not  modeled.  The  times,  temporal  extent,  and  magnitude  of 
travel  time  oscillations  were  deterministically  calculated,  and  are  close  to  the 
measured  values. 

the  PE  to  say  why  those  fluctuations  occurred  and  why  they 
occurred  at  the  times  they  did,  with  the  temporal  extent  and 
amount  of  travel  time  shift  they  have.  Therefore,  ray  tracing, 
although  not  as  reliable,  will  be  used  to  model  the  propagation 
when  NLIWs  are  present.  The  environmental  model  presented 
in  Section  III  will  be  used  to  provide  the  range-  and  time-depen¬ 
dent  environmental  input. 

B.  Ray  Tracing  Results 

When  NLIWs  are  within  the  acoustic  path,  some  of  the  ar¬ 
rivals  show  large  oscillatory  variations,  whereas  the  rest  do  not. 
Ray  tracing  is  used  to  model  the  arrival  time  changes  using  the 
environmental  model  presented  in  Section  III  as  the  range-  and 
time-dependent  environmental  input. 

The  algorithm  used  to  trace  the  rays  is  to  integrate  the 
differential  equations  of  the  rays  with  range  as  the  independent 
variable,  and  with  depth,  vertical  wave  number  divided  by 
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frequency,  and  travel  time  as  the  dependent  variables.  The 
MATLAB  routine  odell3  [11]  performed  the  integration. 

The  sound-speed  field  (Fig.  7)  derived  from  the  CTD  chain 
data  and  interpolated  onto  the  acoustic  path  is  used  as  the  range- 
dependent  sound-speed  field  for  ray  tracing.  The  time  depen¬ 
dence  of  the  sound-speed  field  is  achieved  by  moving  the  en¬ 
tire  sound-speed  field  with  the  wave  speed  0.6447  m/s,  given 
in  Section  III.  Range-dependent  ray  tracing  calculations  were 
carried  out  to  simulate  acoustic  data  between  16:24:00Z  and 
17:26:00Z.  The  arrival  times  for  rays  BT,  TB,  TBT,  and  SBS 
are  simulated.  At  16:44:00Z,  the  source  is  approximately  under 
the  peak  of  the  first  wave. 

The  ray  tracing  results  for  the  four  rays  are  overlaid  on 
acoustic  data  in  Fig.  9.  An  overall  constant  travel  time  of  the 
unmodeled  axial  arrival  used  to  define  the  reduced  travel  time 
was  the  only  adjustable  parameter.  It  is  apparent  that  the  ray 
tracing  results  not  only  successfully  captured  the  times,  extents, 
and  magnitudes  of  the  arrival  time  oscillations  (1  ms)  but  also 
the  gradual  arrival  time  change  of  rays  TBT  and  SBS  due  to  the 
depression  of  the  thermocline.  Two  TBT  rays  existed  most  of 
the  time,  and  the  data  also  often  show  two  separated  arrivals. 

By  closely  examining  the  ray  fields  for  some  specific  cases, 
the  knowledge  of  how  the  rays  were  affected  by  the  two  passing 
waves  was  established  and  a  physical  interpretation  is  provided 
in  the  following  section,  including  why  some  of  the  arrivals 
show  large  oscillations,  but  others  do  not. 

V.  Interpretation 

Both  data  and  modeling  show  that  for  the  TB,  BT,  and  TBT 
arrivals,  fluctuations  by  shortening  the  travel  time  by  1  ms  are 
present,  but  the  SBS  arrival  is  much  less  affected.  The  rays  that 
have  larger  travel  time  fluctuations  are  those  that  turn  (T)  in 
the  strong  thermocline.  The  TBT  rays  have  four  fluctuations  to 
shorter  time,  while  the  TB  and  BT  rays  each  have  two  such 
fluctuations.  All  these  fluctuations  of  rays  with  upper  turning 
points  have  the  shorter  travel  times  when  the  position  of  the 


NLIWs  coincide  with  the  upper  turning  points.  Rays  that  pass 
through  the  thermocline  without  turning,  such  as  the  SBS  rays, 
have  smaller  fluctuations.  We  would  like  to  be  able  to  explain 
the  lowering  of  the  travel  time  by  the  following  argument:  The 
sound  speed  in  the  wave  is  greater  than  it  is  at  the  same  depth 
when  the  wave  is  absent.  Thus,  the  travel  time  along  any  path 
that  passes  through  the  wave  is  shorter  than  the  travel  time  along 
the  same  path  without  the  wave.  We  choose  this  path  pO  to  be  the 
ray  in  the  absence  of  the  wave.  By  Fermat’s  principle,  with  the 
wave  present,  the  actual  ray  pi  has  a  shorter  travel  time  than  the 
path  pO.  The  ray  pi  with  the  wave  has  shorter  travel  time  than 
the  path  pD  with  the  wave,  which  in  turn  has  a  shorter  travel  time 
than  the  ray  pO  without  the  wave.  Therefore,  the  ray  travel  time 
with  the  wave  is  shorter  than  that  without  the  wave.  The  trouble 
with  this  argument  is  that  Fermat’ s  principle  only  requires  the 
ray  to  be  a  (local)  minimum  before  it  reaches  a  caustic.  Beyond 
the  caustic,  it  is  only  a  saddle  point  [12];  there  are  distortions 
of  the  ray  that  give  shorter  travel  times.  Each  of  the  rays  TBT, 
TB,  and  BT  has  at  least  one  caustic.  Thus,  the  argument  must 
be  modified. 

The  following  modified  argument  also  uses  Fermat’s  prin¬ 
ciple,  but  allows  the  ray  travel  time  to  be  a  saddle  point,  rather 
than  only  a  minimum.  When  the  sound- speed  field  changes 
slightly,  the  travel  time  difference  between  the  physical  rays 
in  the  original  and  changed  environments  can  be  evaluated  in 
two  steps.  In  the  previous  argument,  we  first  changed  the  sound 
speed  along  the  original  ray  path,  and  then  changed  the  path 
from  the  original  ray  to  the  ray  for  the  changed  sound- speed 
field.  We  also  consider  the  two  steps  taken  in  the  opposite 
order.  For  small  changes  of  the  sound-speed  field,  the  travel 
time  change  along  a  fixed  path  is  linear  in  the  size  of  the  per¬ 
turbation.  Therefore,  that  step  always  reduces  the  travel  time. 
The  vertical  deviation  Sz(x)  of  a  ray  due  to  a  small  change 
in  the  sound  speed  is  linear  in  the  magnitude  of  that  change, 
because  the  bending  of  the  ray  is  proportional  to  the  gradient 
of  the  sound  speed.  For  either  a  minimum  or  a  saddle  point,  the 
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Fig.  1 1 .  Computed  travel  time  differences  between  rays  when  the  wave  is  absent  and  when  the  wave  is  present,  for  three  ray  species.  Block  “a”  shows  how  to  read 
the  entries  in  blocks  “b,”  “c and  “d.”  Each  travel  time  difference,  shown  on  the  diagonal,  is  decomposed  into  the  sum  of  that  due  to  changing  the  sound-speed 
field  without  changing  the  path,  represented  by  a  horizontal  arrow,  and  that  due  to  changing  the  path  without  changing  the  sound-speed  field,  represented  by  a 
downward  arrow.  Either  of  the  two  partial  changes  can  be  done  first,  leading  to  a  horizontal  arrow  at  the  top  of  the  block  added  to  a  downward  arrow  on  the  right,  or 
a  downward  arrow  on  the  left  added  to  a  horizontal  arrow  on  the  bottom.  Blocks  “b,”  “c,”  and  “d”  are  labeled  by  the  three  ray  species.  Only  one  of  the  two  ordering 
of  the  partial  changes,  which  is  indicated  by  a  dashed  arrow,  has  the  quadratic  path  distortion  effect  much  smaller  than  the  linear  sound-speed  change  effect,  but 
the  other  ordering  has  approximately  equal  partial  changes. 


travel  time  in  a  given  sound-speed  field  is  a  stationary  point  in 
6z(x),  hence  quadratic  in  the  deviation  of  the  path  from  a  ray. 
In  particular,  the  difference  in  travel  times  between  the  rays  for 
two  sound-speed  fields  differing  by  a  small  perturbation,  when 
both  times  are  calculated  with  either  one  of  the  sound-speed 
fields,  is  quadratic  in  the  sound-speed  difference.  For  small 
enough  sound-speed  field  perturbations,  the  linear  effect  of  the 
sound-speed  difference  along  either  curve  is  greater  than  the 
quadratic  effect  due  to  the  two  paths  being  different.  Since  the 
linear  effect  reduces  the  travel  time,  the  total  effect  reduces  the 
travel  time. 

It  only  remains  to  be  seen  whether  the  path  distortion  effect 
is  indeed  small  enough  between  the  absence  and  presence  of 
the  wave,  allowing  one  to  interpret  the  shortened  travel  time 
as  primarily  due  to  the  increased  sound  speed  in  the  wave.  We 
choose  three  cases  to  examine.  There  are  actually  two  TBT  rays. 
Two  of  our  cases  are  the  fluctuations  of  these  two  rays  at  time 
16:48:00Z.  The  third  case  is  the  fluctuation  of  the  TB  ray  at 
16:47:00Z.  These  times  are  those  of  the  minimum  travel  time  in 
the  first  fluctuation,  when  the  first  (or  only)  turning  point  coin¬ 
cides  with  the  first  NLIW.  The  rays  with  and  without  turning  in 
an  NLIW  are  shown  in  Fig.  10. 

The  computed  travel  time  differences  for  these  three  rays  are 
given  in  Fig.  11.  Each  corner  of  each  block  “b,”  “c,”  and  “d” 
represents  the  time  for  a  particular  path  in  a  particular  sound- 
speed  field.  The  upper  corners  are  for  the  path  that  is  the  ray 
when  the  wave  is  absent,  and  the  lower  corners  for  the  path  that 
is  the  ray  with  the  wave.  The  left  corners  are  for  the  sound-speed 
field  without  the  wave,  and  the  right  for  the  sound-speed  field 
with  the  wave.  Thus,  the  upper  left  and  lower  right  corners  are 
the  physical  cases  where  the  paths  are  the  rays.  The  arrows  with 
numbers  indicate  the  time  difference  between  the  corner  at  the 
head  of  the  arrow  and  that  at  the  tail.  The  travel  times  for  all 


the  paths  decrease  (horizontal  arrows)  with  the  wave  present, 
as  they  must,  because  the  sound  speed  everywhere  is  larger  or 
equal  to  that  without  the  wave. 

Let  us  first  examine  the  upper  TBT  ray  in  some  detail.  The 
two  possible  ways  of  taking  the  two  steps  to  change  from  the  ray 
without  the  wave  to  that  with  the  wave  are  shown  in  block  “b” 
of  Fig.  11.  Starting  from  the  top  left  comer  and  moving  verti¬ 
cally  from  downward,  we  distort  the  path.  This  causes  the  travel 
time  to  lengthen  by  0.234  ms.  Now  moving  from  left  to  right, 
the  sound-speed  field  is  changed  to  include  the  waves,  causing 
the  travel  time  to  shorten  by  1 .066  ms,  considerably  greater  than 
the  amount  due  to  path  distortion.  Therefore,  the  modified  argu¬ 
ment  works  here;  the  quadratic  change  due  to  path  distortion 
is  considerably  smaller  than  that  due  to  the  sound-speed  field 
change.  However,  if  we  start  from  the  upper  left  comer,  but 
move  horizontally  first,  the  travel  time  shortens  by  0.460  ms  due 
to  the  change  of  sound-speed  field.  Then,  moving  downward,  the 
travel  time  shortens  further  by  0.373  ms  due  to  path  distortion. 
The  amount  of  change  in  the  two  steps  is  of  the  same  order  of 
magnitude.  The  argument  that  the  quadratic  change  due  to  path 
distortion  is  much  smaller  than  that  due  to  the  sound-speed  field 
change  does  not  work  for  this  order  of  taking  the  two  steps. 

The  lower  TBT  ray  and  the  TB  ray  are  shown  in  blocks  “c” 
and  “d”  of  Fig.  1 1 .  Similar  results  are  obtained  for  these  rays  as 
for  the  upper  TBT  ray,  except  that  the  modified  argument  works 
when  the  sound-speed  change  step  is  taken  before  the  path  dis¬ 
tortion  step.  In  the  order  that  works  for  the  upper  TBT  ray,  the 
two  steps  reduce  the  travel  time  in  roughly  equal  amounts. 

In  each  case,  for  one  of  the  two  orders,  the  change  of  sound 
speed  can  be  considered  small  enough  to  make  the  quadratic 
contribution  small  compared  to  the  linear  contribution.  How¬ 
ever,  for  the  other  order,  the  contributions  are  about  equal;  the 
quadratic  term  is  never  much  greater  than  the  linear  term.  In 
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summary,  the  modified  argument  can  be  taken  in  two  ways:  if 
the  favorable  order  is  taken,  the  argument  works;  otherwise,  the 
argument  is  borderline. 

The  SBS  ray  has  much  smaller  travel  time  fluctuations  (0.3 
ms)  than  the  rays  with  upper  turning  points.  There  are  four 
passes  through  the  thermocline,  and  each  pass  may  hit  one  of 
the  two  waves,  leading  to  eight  opportunities  for  a  travel  time 
decrease.  There  are  five  such  fluctuations  found  in  the  data;  the 
eight  possible  decreases  overlap  in  time  and  the  travel  time  does 
not  return  to  its  no-wave  value.  This  ray  passes  through  the  ther¬ 
mocline  at  an  angle  of  about  12°  and,  because  the  angle  is  large, 
it  encounters  the  wave  for  only  a  short  distance,  so  the  travel 
time  change  is  small. 


VI.  Discussion 

The  nonlinear  internal  waves  observed  with  the  towed  CTD 
chain  measurements  have  been  interpolated  in  time  and  space 
onto  the  acoustic  path  well  enough  to  be  able  to  deterministi¬ 
cally  model  the  travel  time  data  fairly  well.  The  interpolation 
is  not  entirely  straightforward,  as  the  two  sides  of  the  acoustic 
path  were  rather  different.  Assumptions  had  to  be  made,  such  as 
to  what  to  do  about  the  absence  of  the  back  of  the  second  wave 
on  the  southwest  leg  going  seaward,  and  with  the  first  north¬ 
east  leg,  going  shoreward,  not  starting  far  enough  southeast  to 
contain  the  back  of  the  second  wave.  No  attempt  was  made  to 
interpolate  smaller  waves,  which  primarily  influence  the  inten¬ 
sity  rather  than  the  travel  time. 

The  NLIWs  have  a  significant  effect  on  the  travel  time,  even 
though  the  waves  during  this  experiment  were  rather  small  in 
amplitude.  Based  on  experience  from  the  nearby  SWARM  [2] 
experiment,  several  times  larger  waves  were  expected.  If  they 
had  been  bigger,  more  pronounced  effects  would  have  occurred. 
If  the  thermocline  was  depressed  enough  so  that  it  went  below 
the  depth  that  had  been  the  sound  axis,  even  the  direct  path 
would  have  been  strongly  affected.  (This  path  is  a  bunch  of  rays, 
rather  than  one,  at  a  1-km  range  with  the  sound-speed  profile 
that  occurred.) 

Travel  times  decreased  most  notably  when  the  upper  turning 
point  of  a  ray  was  in  (or  near)  the  thermocline  and  an  NLIW 
encountered  this  upper  turning  point.  Every  time  this  encounter 
occurred,  the  travel  time  decreased  by  about  1  ms.  The  presence 
of  caustics  invalidates  a  straightforward  Fermat’s  principle  in¬ 
terpretation  with  a  minimum  travel  time  for  the  ray.  Our  attempt 
to  interpret  this  systematic  behavior  based  on  Fermat’s  principle 
of  stationary,  rather  than  minimum  travel  time,  where  ray  path 
change  effects  are  assumed  small  was  only  partly  successful. 
A  two-step  process  was  used  to  move  from  the  ray  without  the 
wave  to  that  with  the  wave.  The  effect  of  ray  path  change  de¬ 
pended  on  the  order  of  the  two  steps.  Only  one  order,  and  not 
always  the  same  one,  resulted  in  a  small  travel  time  difference 
due  to  the  distortion  of  the  ray  path  by  the  NFIW.  For  the  other 
order,  it  was  comparable  to  the  travel  time  decrease  due  to  the 
higher  sound  speed  along  the  ray  path. 

The  study  has  only  been  concerned  with  the  direct  problem. 
The  ability  to  carry  out  the  inverse  problem  would  depend  on 


what  is  assumed  known.  With  only  the  acoustic  data,  it  is  un¬ 
likely  that  a  meaningful  inverse  could  be  done.  If  the  acoustics 
and  the  CTD  profile  from  15:30:00Z  were  available,  there  would 
be  a  much  better  chance  of  useful  inversion,  although  the  deep¬ 
ening  of  the  thermocline  over  the  hour  from  the  time  of  the  CTD 
would  have  to  be  allowed  for  in  the  inversion  method. 
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Abstract:  Downward  looking  sonar,  such  as  the  chirp  sonar,  is  widely 
used  as  a  sediment  survey  tool  in  shallow  water  environments.  Inversion 
of  geo-acoustic  parameters  from  such  sonar  data  precedes  the  availabil¬ 
ity  of  forward  models.  An  exact  numerical  model  is  developed  to  initiate 
the  simulation  of  the  acoustic  field  produced  by  such  a  sonar  in  the  pres¬ 
ence  of  multiple  rough  interfaces.  The  sediment  layers  are  assumed  to 
be  fluid  layers  with  non-intercepting  rough  interfaces. 
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1.  Introduction 

Chirp  sonar  is  widely  used  as  a  survey  tool  to  image  shallow  water  sediment  structure. 
Usually,  it  has  a  beam  pattern  aimed  vertically  toward  the  seafloor  and  operates  at  a 
nominal  center  frequency  of  several  kilohertz  with  a  wide  frequency  band  on  the  order 
of  10  kHz,  providing  a  depth  resolution  on  the  order  of  10  cm.  Chirp  sonar  data  inver¬ 
sion  has  been  used  to  determine  sediment  layering  and  physical  properties  in  underwater 
acoustics  applications.1-4  A  curious  fact  is  that  an  exact  forward  model  predicting  the 
received  sound  field  by  such  a  sonar  for  a  given  physical  sediment  structure  has  not 
been  developed.  The  availability  of  such  predictive  models  would  help  the  interpretation 
of  chirp  sonar  data  and  validation  of  geo-acoustic  inversion  results  based  on  such  data. 
This  paper  is  an  initial  step  toward  such  a  forward  model  and  strives  for  an  exact  solu¬ 
tion  of  interface  roughness  scattering  using  the  integral  equation  method.  It  is  not  meant 
for  direct  application  to  real  sonar  systems,  but  to  provide  a  basis  on  which  practical 
models,  such  as  the  recent  model  based  on  a  transfer  function  approach,5  can  be 
validated. 

The  sediment  is  modeled  in  the  two-dimensional  space  and  consists  of  a  set  of 
fluid  layers  separated  by  one-dimensional  rough  interfaces.  Each  layer  is  assumed  to  be 
homogeneous  with  constant  sound  speed,  density,  and  attenuation  coefficient.  The 
rough  interfaces  are  assumed  to  be  single-valued  and  never  to  cross  one  another,  but 
otherwise  allowed  to  be  arbitrary.  This  initial  model  does  not  allow  for  sediment 
volume  heterogeneity.  In  addition  to  sediment  acoustics,  this  model  has  potential  appli¬ 
cation  in  imaging  Arctic  ice  from  underneath  the  ice  layer.  When  the  roughness  is 
small  compared  to  the  acoustic  wavelength,  perturbation  method  can  be  used  to  treat 
multiple  rough  layers.6  Here  we  do  not  make  such  an  assumption  and  the  rough  sur¬ 
face  amplitudes  are  allowed  to  be  large  as  compared  to  the  acoustics  wavelength. 

2.  Single  rough  interface  problem 

Let  the  one-dimensional  rough  interface  £(x)  separate  two  fluid  half-spaces.  The  coor¬ 
dinate  system  is  defined  such  that  the  positive  z-axis  points  upward  and  the  x-axis  to 
the  right.  A  monochromatic  incident  wave  /?inc(r)  with  frequency  /  impinges  from 
above  the  interface,  where  r  =  (x,z).  The  upper  medium  has  sound  speed  c i,  density  pi, 
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and  attenuation  factor  S\,  which  is  the  ratio  of  the  imaginary  part  to  the  real  part  of 
the  wavenumber  k\  =  (2nf/ci)(\  +i5 i),  and  the  lower  medium  has  corresponding  sound 
speed  c2,  density  p2,  and  attenuation  factor  d2,  as  well  as  wavenumber  k2.  Following 
Thorsos,7,8  the  coupled  integral  equations  with  two  unknown  quantities  p\(r)  and 
on  the  rough  interface  are 


Pi  (r)  =  2/?„ 


*«  =  2 


;(r)  -  2 


Ho(ki\r  —  r'l) 


P2^Pl(0 

Pi  dn' 


dpi  (r7)  dHo(ki\r  —  r'|) 
dn '  dn' 

dHl0(k2  |r  — r'|) 


Pi(r') 


ds' 


dn' 


p  i(0 


(i) 


ds' , 


where  d/dn'  is  the  surface  normal  derivative,  and  Hq  is  the  Hankel  function  of  the 
first  kind.  The  integrals  are  along  the  rough  interface  and  all  quantities  in  Eq.  (1)  are 
evaluated  on  the  interface.  The  unknowns  in  Eq.  (1)  are  chosen  to  be  the  field  and  its 
normal  derivative  in  the  upper  medium.  Corresponding  quantities  in  the  lower  medium 
are  found  through  the  boundary  conditions 


[Pi(r)  =Pi{  r) 

|  1  <9pi(r)  _  1  dp2(r)  (2) 

[  Pi  dn  p2  dn 


When  pi( r)  and  are  solved,  the  scattered  field  anywhere  in  media  1  and  2  are 

given  by 


Pjr)=Pi(r)-Pinc(r)=4 


|r-r'|) 


,  api(r')  dH^\kx  |r-r'|) 


dn' 


pi(’)='~A\ 


a"“'l(Sr  ~  -P«) 


Pi  dn' 


dn' 


dn' 

ds’. 


p\  (O 


ds' 


(3) 


By  discretizing  the  interface  into  N  equal  segments  Ax  =  LIN,  where  L  is  the  total 
length  of  the  ensonified  interface,  the  integral  Eq.  (1)  is  recast  into  the  matrix  form 

q  =  Ay,  (4) 


where  q  is  a  known  column  vector  of  length  IN.  The  y'th  element  of  q  is 


„  =  /pinc(r/),  j=\,2,...,N 
%  1 0,  j=(N+l),(N  +  2),...,2N, 


(5) 


and  y  is  an  unknown  column  vector  of  length  2 TV.  The  rath  element  of  y  is 


ym  — 


Axym  dpi  (rm') 


4  i 

A-*7m 

4  / 


dn' 

P\  (jm  )? 


ra  =  1,2,...,  N 

ra=(7V+l),(7V  +  2),...,27V, 


(6) 


where  ym  =  y  1  +  ^q.  (4),  A  is  a  2 N  x  2 N  matrix,  which  can  be  more  con¬ 
veniently  expressed  as  four  NxN  sub-matrices:  A=  (^2!  ^22),  where  the  element  at 

the  jth  row  and  rath  column  of  each  of  the  sub-matrices  is,  respectively, 
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H^(kx\rj-rm'\),  j=l,2,...,N-  m  =  1,2,..., TV;  (j  ^  m) 


An  - 
Ajm  ~ 


A21  -  l 

jm  ~  ' 


H^\kxAxyj/{2e)),  j  =  1,2,..., TV;  (j  =  m ), 

'  ^H^(k1\rj-rmt\)  j  =  1,2,..., JV;  m  =  1,2 . N-  (j^m) 
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Pi 


r  dH^ [){ki\rj-Tmf\) 
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,  j  =  1,2,..., m  =  1,2,...,  TV;  (/  +  m)  (7) 


y  =  u,...,W;  y,m), 

Try)  J 


421 

jm 


dH^ 


j=  m=  1,2,..., N;  (j  ±  m) 


Y^r+(^\ 

\yjAx  ny]  1 


The  matrix  Eq.  (4)  can  be  solved 


y  =  Sq,  (8) 

where  S  =  A  1 .  The  matrix  A,  hence  S,  is  a  function  of  the  rough  interface  and  media 
properties  above  and  below  the  interface,  and  is  independent  of  the  incident  wave 
Pi nc(r),  or  q.  This  property  will  be  exploited  in  problems  of  scattering  from  multiple 
rough  interfaces.  Another  useful  property  of  the  matrix  A  is  how  the  matrix  Eq.  (4) 
would  change  when  the  incident  wave  impinges  from  below  the  interface.  In  Eq.  (1), 
interchange  the  subscripts  1  and  2,  and  still  express  pi( r')  and  dp^J  ^  as  the  unknowns 

using  boundary  conditions  (2),  one  finds  that  the  matrix  equation  for  an  incident  wave 

coming  from  the  lower  medium  has  the  same  matrix  form  as  Eq.  (4) 

q  =  Ay,  (9) 


where  the  elements  of  q  are  similar  to  q,  but  its  second  half  is  the  incident  wave 


%■ 


0, 

Pine  (r/)  5 


jm  1,2,...,  A 

j={N+  1),(A  +  2),...,2A. 


(10) 


This  symmetry  property  means  that  the  scattering  matrix  S  uniquely  governs  sound 
scattering  from  a  particular  rough  interface  for  incident  waves  coming  either  from 
above  or  below.  Once  S  is  obtained,  scattering  from  the  interface  for  a  given  frequency 
for  any  normal  incident  wave  is  uniquely  determined. 

Once  y  is  known,  the  scattered  field  anywhere  in  space  can  be  calculated  using 
Eq.  (3),  or  its  equivalent  matrix  form.  For  M  observation  points  in  the  upper  medium, 
the  scattered  field  is 


qi  =  Twy, 


(ID 
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where  qi  is  a  vector  of  length  2 M,  whose  first  M  elements  are  zero,  and  the  second  M 
elements  are  the  scattered  field  /?i5(rm)  evaluated  at  the  M  points  rm,  m=  1,2 in 
medium  1 .  The  upward  propagation  matrix  Tw  is  a  2 M  x  IN  matrix  whose  first  M 
rows  are  zero,  and  the  rest  are 


1  umj 


-dH*)^m~rj  ^,m  =  (M+l),(M+2),...,2M,j=(N+l),(N  +  2),...,2N. 

(12) 


Similarly,  for  any  M  observation  points  in  the  lower  medium,  the  scattered  field  is 


q2  =  T^y, 


(13) 


where  q2  is  a  vector  of  length  2 M,  whose  second  M  elements  are  zeros,  and  whose  first 
M  elements  are  the  scattered  field  /?2(rm)  evaluated  at  the  M  points  rm,  m  =  1,2 in 
medium  2.  The  downward  propagation  matrix  Td  is  also  a  2 Mx  2 N  matrix,  though  its 
second  M  rows  are  zero,  and  the  first  M  rows  are 


Tdmi  — 


-^H{0l\k2 |rm-r/|),  m  =  1,2,..., M,  j=l,2,...,N 
Pi 


dH^  (k2\rm  -  tj'\) 


(14) 


dn! 


m  =  1,2 j  =  (TV  +  1),  (7V  +  2),  ...,27V. 


3.  Multiple  rough  interfaces  problem 

For  a  single  rough  interface  we  have  shown  that  each  rough  interface  £(x)  has  an  asso¬ 
ciated  scattering  matrix  S  that  is  a  function  of  the  parameters  of  the  media  above  and 
below  the  interface,  and  uniquely  determines  the  scattered  field  for  any  incident  field. 
Two  propagation  matrices,  Tu  and  Td,  allow  calculation  of  the  scattered  field  above  or 
below  the  rough  interface  once  the  field  and  its  normal  derivative  on  the  interface  are 
calculated.  When  there  are  two  or  more  rough  interfaces,  the  previous  results  are  used 
to  construct  the  scattered  field. 

In  a  scenario  of  two  rough  interfaces,  let  (i(v)  be  a  rough  interface  situated 
above  another  rough  interface  £2(x).  Also  assume  that  the  two  interfaces  never  intersect 
with  each  other,  or  Ci(v)  >  (2(X)-  While  the  parameters  above  and  below  Ci(-V)  are  the 
same  as  those  for  a  single  rough  interface,  the  medium  beneath  C2(v)  has  sound  speed 
c3,  density  p3,  attenuation  factor  <53,  as  well  as  wavenumber  k3.  Each  rough  interface 
has  an  associated  scattering  matrix  and  a  pair  of  propagation  matrices,  and  are  given 
the  names  S y,  and  Tju  and  y=l,2.  The  discretized  unknown  field  and  its  normal 
derivative  for  each  rough  interface  are  always  chosen  to  be  the  ones  above  the  corre¬ 
sponding  interface.  They  are  designated  as  y3  and  y2  and  are  expressed  in  the  same 
manner  as  those  given  in  Eq.  (6)  for  each  interface.  For  interface  Ci(x),  there  are  two 
waves  impinging  onto  it:  one  is  the  original  incident  wave  q  coming  from  above,  the 
other  is  T2wy2  that  is  scattered  upward  by  Ciix)  from  below.  Whereas  for  C2(X)5  there  is 
only  one  incident  wave,  T^yi,  forward  scattered  by  Ci(x).  Using  Eqs.  (8)  and  (9)  and 
linear  superposition  yields 

fyi  =S1q  +  S1T2„y2  (J5) 

ty2  =  SiTidyj. 

Solving  the  equations  for  the  unknown  on  the  first  interface  yields  y!  =  Siq 
+  (SiT2MS2Tlflj)y1?  or 
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yi=T — Siq,  (16) 

1  -  -*M,2 

where  Rij2  =  SiT2mS2T1^  is  the  multiple  scattering  matrix  for  the  two-layer  system.  The 
result  in  Eq.  (16)  has  a  “ray”  interpretation  and  it  can  be  more  explicitly  seen  if  the 
matrix  in  the  denominator  is  expressed  using  the  binomial  expansion 

i  oo 

yi=r-S-S1q  =  ^R^2S1q,  (17) 


where  the  n  =  0  term  is  for  the  case  where  the  second  interface  is  absent,  and  the  n  =  1 
term  is  used  where  the  second  interface  scatters  once.  In  general,  the  nth  term  is  for 
the  case  where  the  wave  has  been  reflected  between  the  two  interfaces  n  times. 

For  N  rough  interfaces,  a  procedure  similar  to  Eq.  (15)  can  be  used  to  con¬ 
struct  the  multiple  scattered  field.  In  particular,  if  one  is  interested  in  finding  the  field 
scattered  back  to  medium  1,  as  would  be  used  for  chirp  sonar,  a  general  multiple  scat¬ 
tering  matrix  R1>N  can  be  derived  to  be 


yi  = 


i 

1  —  Ri,v 


Siq, 


(18) 


where  the  multiple  scattering  matrix  is  given  recursively 


Rl  ,N  = 


Ri,v-i 
1  -  Rjv-i,# 


(19) 


For  example,  when  there  are  three  rough  interfaces 

^  _  Rl2  _  SiT2mS2Ti</ 

1,3  “  1  -  R2l3  “  1  -  s2t3„s3t2/ 


(20) 


4.  Numerical  example 

To  demonstrate  the  application  for  two  rough  interfaces,  consider  a  sound  source  and 
receiver  co-located  20  m  above  the  mean  seafloor.  The  point  source  emits  1  ^iPa  at  1  m 
range  and  has  a  Gaussian  beam  pattern  such  that  the  incident  sound  amplitude  on 
the  mean  seafloor  drops  to  l/e  at  5.6m  from  the  center  of  the  beam.  The  source  has  a 
center  frequency  of  3500  Hz,  and  a  Gaussian  spectrum  such  that  the  baseband  has 
half  power  at  1000  Hz.  The  sediment  is  a  layer  of  sand  with  a  mean  thickness  of  3  m, 
which  is  in  turn  on  top  of  a  mud  half  space.  The  parameters  for  the  three  layers  are 
given  in  Table  1.  Both  interfaces  have  the  same  power  spectrum  as  provided  by 

the  ONR  Reverberation  Modeling  Workshop:9,10  P  =  n^2^K2 y  with  h  =  0.316m, 

Kl  =  2.510“ 3  m_1  for  the  water-sand  interface  specified  as  “typical,”  and  A  =  0.1m, 
Kl=  lm_1  for  the  sand-mud  interface  specified  as  “rough.”  The  scattered  fields  are 
calculated  at  a  set  of  172  frequencies  around  the  center  frequency  with  37.5-Hz  incre¬ 
ments;  Fourier  synthesis  is  used  to  obtain  time  domain  fields.  In  this  example  (Fig.  1), 


Table  1.  Medium  parameters. 


C  (m/s) 

p  (kg/m3) 

S 

Water 

1500 

1000 

0 

Sand 

1740 

2000 

0.01 

Mud 

1485 

1400 

0.001 
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Fig.  1.  (Color  online)  The  first  three  arrivals:  the  first  is  from  the  water-sand  interface  with  roughness  described 
by  the  “typical”  spectrum,  the  second  from  the  sand-mud  interface  with  roughness  described  by  the  “rough” 
spectrum,  and  the  third  also  from  the  sand-mud  interface  with  a  round  trip  between  the  two  interfaces.  The 
dashed  curves  are  for  flat  interfaces.  The  thin  solid  curves  are  arrivals  from  seven  rough  interfaces. 

there  are  three  arrivals.  The  first  is  the  scattered  field  by  the  water-sand  interface 
(blue),  the  second  is  (a)  forward  scattered  by  the  water-sand  interface  onto  the  sand- 
mud  interface,  (b)  backscattered  to  the  water-sand  interface,  and  then  (c)  forward  scat¬ 
tered  to  the  water  column  (green).  The  third  is  similar  to  the  second  but  with  an  addi¬ 
tional  round  trip  inside  the  sand  layer  (red).  Results  for  the  flat  interfaces  are  also 
given  (Fig.  1,  dashed)  which  has  a  16dB  loss  from  propagation  from  the  20  m  height. 
The  main  effects  of  the  roughness  as  compared  to  the  flat  interface  are  (1)  to  cause  the 
main  return  to  fluctuate  as  a  function  of  position,  and  (2)  to  produce  a  coda,  or  scat¬ 
tered  wave  that  follows  the  main  arrival.  Later  arrivals  have  more  fluctuation  due  to 
multiple  scattering.  Figure  2  shows  results  over  a  range  of  100  m  containing  100  pings, 
from  which  several  observations  can  be  made:  (i)  the  low  spatial  frequency  roughness 
of  the  interfaces  can  be  traced  by  the  scattered  intensities;  (ii)  the  coda  appears  as 
“noise”  in  between  the  layers,  not  unlike  those  observed  in  chirp  data;  and  (iii) 
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Fig.  2.  Simulation  for  two  rough  interfaces.  There  are  100  pings  over  a  200  m  range.  The  mean  distance 
between  the  two  interfaces  is  3  m.  The  colorbar  shows  intensity  in  dB.  The  faint  signals  arriving  earlier  than  the 
water-bottom  interface  (time  earlier  than  26  ms)  are  numerical  artifacts  due  to  finite  bandwidth  in  the  Fourier 
synthesis. 
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multiple  scattered  arrivals  are  appreciable,  though  about  20  dB  lower  than  those  scat¬ 
tered  by  the  lower  interface. 

5.  Summary 

This  paper  provides  a  method  for  calculating  exact  scattered  fields  from  a  system  of 
rough  interfaces  in  two-dimensional  space  for  normal  incident  sonar.  It  is  intended  as 
a  basis  of  comparison  for  approximate  methods  such  as  the  Kirchhoff  approximation. 
It  is  also  a  tool  for  studying  statistics  of  scattered  waves  and  validating  techniques  for 
estimating  sediment  properties  using  such  sonar. 
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A  rough-interface  reverberation  model  is  developed  for  range-dependent  environments.  First-order 
perturbation  theory  is  employed,  and  the  unperturbed  background  medium  can  be  layered  and  het¬ 
erogeneous  with  arbitrary  range  dependence.  To  calculate  the  reverberation  field,  two-way  forward 
scatter  due  to  the  slowly  changing  unperturbed  environment  is  handled  by  fast  numerical  methods. 
Backscatter  due  to  small  roughness  superimposed  on  any  of  the  slowly  varying  interfaces  is  handled 
efficiently  using  a  Monte  Carlo  approach.  Numerical  examples  are  presented  to  demonstrate 
the  application  of  the  model.  The  primary  purpose  of  the  model  is  to  incorporate  relevant  physics 
while  improving  computational  speed.  ©  2012  Acoustical  Society  of  America. 
[http://dx.doi.Org/10.1121/l.4707437] 
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I.  INTRODUCTION 

This  paper  expands  model  capabilities  for  range- 
dependent  shallow  water  reverberation  due  to  roughness  of 
interfaces,  including  the  sea  surface  and  water-seafloor 
interface  as  well  as  buried  interfaces  such  as  the  sediment 
basement.  In  such  problems,  scattering  is  diffusely  distrib¬ 
uted  in  space  and  contributes  to  background  reverberation 
as  compared  to  distinct  returns  due  to  targets  and  target- 
like  clutter.  The  practical  issue  challenging  such  modeling 
efforts  is  fidelity  vs  speed.  When  the  environment  is  speci¬ 
fied,  the  reverberation  problem  can  in  principle  be  solved 
exactly  when  sufficient  computational  power  is  available 
using  numerical  methods,  such  as  the  finite  element 
method.  However,  such  brute  force  approaches  cannot 
presently  render  practically  useful  solutions.  To  obtain 
practical  solutions,  compromise  is  needed  between  compu¬ 
tation  speed  and  accuracy.  In  the  foreseeable  future,  such  a 
compromise  will  be  the  primary  modus  operandi.  Then  the 
questions  are:  First,  with  given  mission-required  fidelity 
versus  uncertainty,  is  there  an  approximate  model  that  can 
predict  reverberation  with  sufficient  speed?  Second,  what 
is  the  uncertainty  in  the  reverberation  prediction?  The 
primary  goal  of  this  paper  is  to  address  the  first  question. 
In  particular,  this  article  develops  a  physics-based  model 
that  calculates  reverberation  in  a  waveguide  with  arbitrary 
range-dependent  bathymetry  and  stratification  where  the 
interfaces  have  randomly  distributed  roughness.  The 
small-roughness  perturbation  approximation  is  used  in 
which  the  roughness  responsible  for  acoustic  backscatter- 
ing  is  assumed  to  be  much  smaller  than  the  acoustic  wave¬ 
length.  It  is  important  to  note  that  the  smallness  criterion 
does  not  apply  to  larger-scale  bathymetry.  On  the  one 
hand,  the  range-dependent  propagation  is  handled  by  fast 
numerical  methods  such  as  the  parabolic  equation  (PE), 
hence  the  computational  speed  is  reasonably  fast;  on  the 
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other  hand,  because  the  basis  of  this  approach  is  physical 
rather  than  phenomenological,  the  results  can  be  used  to 
assess  the  appropriateness  of  other  faster  but  possibly  less 
accurate  methods. 

Except  for  exact  numerical  approaches  such  as  the  finite 
element  method,1  which  have  been  made  available  with 
increased  computer  capacity,  reverberation  models  assume 
the  reverberant  return  is  a  result  of  two-way  outgoing-and 
return-propagation  with  a  single  backscatter  coupling  the 
two.  The  outgoing-and  return-propagation  may  itself  involve 
a  different  kind  of  scattering  in  the  sense  that  a  range- 
dependent  environment  alters  the  waveguide  field  compared 
to  the  range-independent  case,  but  scattering  that  reverses 
the  direction  of  propagation  is  assumed  to  be  negligible.  In 
the  following,  the  term  forward  scattering  is  applied  to  com¬ 
putation  of  both  the  outgoing  and  return  field.  This  multiple 
forward  scatter,  single  backscatter  assumption  is  reasonable 
except  in  cases  where  the  mechanism  responsible  for  back¬ 
scatter  also  causes  an  appreciable  change  in  forward  scatter, 
e.g.,  a  dense  bubble  cloud. 

In  many  cases,  outgoing-and  return-propagation  is  fur¬ 
ther  approximated  by  range -independent  propagation  with 
forward  scatter  ignored.  One  of  the  earliest  and  most  influen¬ 
tial  models  is  by  Bucker  and  Morris,2  who  suggested  using 
normal  mode  intensity  to  approximate  the  two-way  out¬ 
going-and  return-propagation.  In  their  approach,  the  single 
bottom  backscatter,  which  is  specified  by  a  scattering  cross 
section  with  Lambert’s  law  as  its  default  form,  causes  cou¬ 
pling  between  outgoing  and  return  modes  with  incident  and 
scattered  angles  of  the  cross  section  related  to  the  mode  hori¬ 
zontal  wavenumber.  They  further  assumed  that  cross-mode 
contributions  are  unimportant.  When  the  mode  group  speeds 
are  approximated  by  the  reference  water  sound  speed,  the 
computational  speed  is  fast.  This  model  has  been  improved 
in  its  fidelity  in  accuracy,  with  modest  cost  in  speed,  to  allow 
the  individual  normal  modes  to  travel  with  their  own  group 
speeds.3-6 

Specific  scattering  mechanisms,  such  as  bottom 
roughness7-10  and  sediment  volume  heterogeneity11  have 
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been  incorporated  into  the  model  discussed  in  the  preceding 
text,  further  improving  fidelity.  In  these  improvements,  the 
propagation  is  calculated  by  normal  mode  codes  such  as 
KRAKEN12  or  wavenumber  integration  codes  such  as 
OASES.13  In  the  same  category  and  widely  used,  ray  tracing 
is  sometimes  employed  to  calculate  the  propagation.14  This 
offers  competitive  speed,  but  issues  encountered  in  using 
rays  have  to  be  handled  carefully. 

An  especially  fast  method  is  based  on  an  energy  flux 
approximation.3,15-18  The  high  speed  of  this  approach  origi¬ 
nates  from  its  closed-form  nature  when  certain  assumptions 
are  made  regarding  the  boundary  reflection  coefficient,  water 
column  sound  speed,  and  scattering  cross  sections.  This 
method  may  remain  the  preferred  choice  for  applications 
when  speed  is  paramount. 

In  modeling  reverberation  due  to  interface  roughness, 
the  perturbation  method  is  effective  if  the  roughness  devia¬ 
tion  from  a  mean  interface  is  small  compared  to  the  acoustic 
wavelength  and  if  the  maximum  slope  of  the  interface  is  suf¬ 
ficiently  small.19,20  The  roughness  is  commonly  quantified 
by  a  power  spectrum,  and  the  rough  interface  is  assumed  to 
be  spatially  stationary.  In  reverberation  models  using  the 
perturbation  approximation,  the  mean  interface  is  typically 
assumed  to  be  flat;  low  spatial  frequency  variations  do  not 
exist.  In  terms  of  the  roughness  power  spectrum,  the  assump¬ 
tion  is  that  there  is  a  low-wavenumber  scale  below  which  the 
power  spectrum  is  zero.  The  interface  roughness  spectrum 
is  likely  to  be  an  approximate  power-law  over  scales  from 
1  cm  to  1  m  with  no  clear-cut  division  into  low-and  high- 
wavenumber  parts21  (Chap.  6).  Measurements  extending 
over  the  scales  of  interest  in  low-  and  medium-frequency 
reverberation  have  yet  to  be  made,  so  the  question  of  separa¬ 
tion  of  large  and  small  scales  remains  open  in  these  cases. 

In  this  paper,  it  is  assumed  that  the  interface  roughness 
has  two  separable  scales,  one  low- wavenumber  component 
that  may  correspond  depth  changes  that  are  much  larger  than 
the  acoustic  wavelength  but  that  are  slowly  undulating  in 
space;  the  other  high-wavenumber  component  corresponds 
to  roughness  having  vertical  scale  that  is  small  compared  to 
the  acoustic  wavelength  and  horizontal  scale  comparable  to 
the  wavelength.  This  small-scale  roughness  is  superimposed 
on  the  slowly  changing  low-wavenumber  background.  One 
example  is  small  roughness  on  a  sloping  bottom.  With  the 
two-scale  assumption,  it  is  further  assumed  that  the  slowly 
changing  component  causes  forward  scatter  but  no  backscat- 
ter.  The  single  backscatter  is  handled  using  the  first-order 
perturbation  approximation  against  this  slowly  changing 
background.  With  these  assumptions,  the  forward  scatter  in  a 
gently  changing  environment  can  be  handled  by  choosing 
any  efficient,  one-way  propagation  method,  such  as  PE, 
coupled  modes,22  or  the  virtual  source  method,23  making 
modeling  reverberation  in  range-dependent  environments 
practical  even  when  mode  coupling  is  significant.  However, 
the  validity  of  the  two-scale  assumption  is  not  in  general  jus¬ 
tified,  and  the  user  should  be  careful  in  applying  the  method 
to  new  scenarios. 

The  formulation  to  be  given  here  of  first-order  perturbation 
theory  with  a  general,  non-flat  background,  is,  to  the  authors’ 
knowledge,  a  new  result.  In  contrast  to  most  published  work, 


this  formulation  offers  first-order  perturbation  results  on  indi¬ 
vidual  surface  realizations,  rather  than  the  ensemble-averaged 
intensity. 

The  structure  of  the  paper  is  such  that  the  full  acoustics 
problem,  including  the  definition  and  the  main  results,  is  pre¬ 
sented  in  Sec.  II,  while  the  details  of  the  derivation  are  given 
in  Appendix  A.  Section  III  gives  a  range-independent  exam¬ 
ple  for  a  two  half-space  problem  where  the  numerical  result 
for  the  scattering  cross  section  is  compared  to  results  using 
known  theoretical  formulas.  Detailed  steps  of  computation, 
which  are  common  to  other  sections,  are  provided.  In  many 
applications,  the  phenomenological  Lambert’s  law  is  used  as 
the  scattering  cross  section.  In  Sec.  IV,  the  new  formula  is 
applied  to  a  layered  environment  that  has  been  tailored 
to  mimic  Lambert’s  law  over  limited  angular  and  frequency 
ranges.  This  offers  a  step  toward  comparison  of  results 
between  fast,  application  driven  codes  and  physics-based 
methods.  Section  V  is  devoted  to  two  reverberation  problems 
in  which  the  reverberation  time  series  is  of  primary  interest. 
The  first  example  is  reverberation  in  a  Pekeris  waveguide, 
where  other  known  numerical  results  are  available  for  com¬ 
parison.  The  second  example  is  reverberation  in  a  wedge- 
shaped  waveguide,  demonstrating  the  new  formulation  in  a 
range-dependent  environment. 

II.  FORMAL  RESULTS 

The  primary  goal  of  this  work  is  to  develop  a  numeri¬ 
cally  efficient  model  for  reverberation  in  range-dependent 
waveguides,  applicable  at  low  to  mid-frequencies  (nominally 
200  Hz  to  10  kHz).  Attention  will  be  restricted  to  scattering 
by  roughness  of  the  water- seafloor  interface  as  well  as  by 
buried  rough  interfaces,  such  as  the  basement.  Both  the  sedi¬ 
ment  and  basement  (if  any)  will  be  treated  as  acoustic  fluids, 
i.e.,  shear  effects  will  not  be  included  in  the  model.  A  two- 
scale  approach24-26  will  be  used  in  which  the  relief  of  the 
seafloor  is  expressed  as  the  sum  of  a  slowly  varying,  large- 
scale  part  and  a  small-scale  part.  It  is  assumed  that  the 
slowly  varying  part  of  the  seafloor  relief  (bathymetry)  is 
known  and  that  a  propagation  solution  is  available  via  some 
efficient  method  such  as  the  parabolic  equation.  This  will  be 
referred  to  as  the  “zeroth-order”  solution.  This  part  of  the 
problem  incorporates  range  dependence,  and  the  variation  of 
water  depth  and  seafloor  properties  may  be  as  extreme  as  the 
numerical  solution  method  allows.  The  small-scale  part  is  re¬ 
sponsible  for  acoustic  scattering  and  will  be  treated  using 
first-order  perturbation  theory.  This  will  be  referred  to  as  the 
“first-order”  solution. 

The  primary  difference  between  the  present  model  and 
previous  work  lies  in  the  manner  in  which  the  large-scale 
relief  is  modeled  and  in  the  manner  in  which  the  large  and 
small  scales  interact.  For  example,  the  parabolic  equation 
method  has  been  applied  to  calculate  reverberation  due  to 
bottom  roughness,27  where  the  details  of  the  roughness  are 
taken  into  account  by  taking  very  small  steps,  hence  this 
approach  is  not  numerically  efficient.  In  most  applications  of 
the  two- scale  (or  composite-roughness)  approximation,  both 
large-  and  small-scale  relief  are  taken  to  be  random  with  the 
the  effect  of  large  scales  treated  statistically.  The  interaction 
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of  the  large  and  small  scales  is  handled  in  an  ad  hoc  fashion 
with  the  large  scales  assumed  to  simply  tilt  the  landscape 
and,  in  some  cases,  cause  shadowing.  In  the  method  of  this 
article,  a  two-scale  approach  is  developed  from  first  princi¬ 
ples,  and  no  ad  hoc  rules  are  invoked.  The  price  paid  for  this 
more  physically  based  approach  is  that  formal-average 
results  are  not  obtainable,  and  Monte  Carlo  methods  must  be 
used  with  the  small-scale  relief  drawn  from  a  random 
ensemble. 

The  method  of  this  article  faces  the  same  issue  of  sepa¬ 
ration  of  scales  that  has  confronted  two-scale  models  in  the 
past.  Scale  separation  can  be  viewed  as  a  filtering  operation 
in  which  a  low-pass  filter  yields  the  large-scale  relief,  and  a 
high-pass  filter  (using  the  same  cutoff)  yields  the  small-scale 
relief.  For  a  given  scenario,  one  must  find  a  cutoff  such  that 
both  the  zeroth-order  and  first-order  solutions  are  accurate. 
In  some  problems,  no  such  cutoff  exists,  but  when  it  does 
exist,  the  proper  choice  of  cutoff  frequency  becomes  the  pri¬ 
mary  issue.  Once  this  choice  is  made,  the  final  results  should 
be  somewhat  insensitive  the  exact  value  employed.  The 
authors  are  unaware  of  any  set  of  simple,  general  criteria 
that  will  guide  the  choice  of  cutoff  for  the  various  applica¬ 
tions  envisioned.  The  scale  separation  problem  is  not  at  issue 
in  this  article  as  the  zeroth-order  interface  for  all  problems 
considered  is  flat,  although  it  is  sloping  in  one  example. 

The  general  problem  is  illustrated  in  Fig.  1.  A  point 
source  is  situated  at  position  r0,  and  one  wishes  to  predict 
the  bistatic,  average  reverberation  intensity  time  series  at 
position  r.  This  section  will  assume  a  single-frequency 
source,  and  the  issues  of  Fourier  synthesis  of  time  series  and 
Monte  Carlo  averaging  will  be  deferred  to  later  sections.  The 
source  is  assumed  to  produce  unit  pressure  at  unit  distance, 
and  the  pressure  at  the  general  point  r  is  the  Green’s  func¬ 
tion,  satisfying  the  following  equation: 


P(  r)V 


1 

_P(r) 


VG(r,  r0) 


=  — 471(5 (r  -  r0). 


+  k2(r)G(r,  r0) 


(1) 


FIG.  1.  A  range-dependent  waveguide  with  small-scale  perturbation  in 
bathymetry.  The  solid  curve  indicates  the  smooth  part  of  the  bathymetry  for 
which  a  field  soution  is  presumed  known  in  the  form  of  the  Green’s  function 
G(0)(r?  ro)-  The  shaded  regions  define  the  “perturbed”  problem  with  the 
interface  between  dark  and  light  shades  representing  the  “true”  bathymetry. 
Although  not  shown  in  the  figure,  both  the  water  column  and  the  seafloor 
may  be  heterogeneous  as  long  as  these  heterogeneities  are  included  in  the 
calculation  of  G(0)( r,  r0). 


The  density,  p(r),  and  wavenumber,  £(r),  can  be  position  de¬ 
pendent,  and  attenuation  is  included  by  allowing  k(r)  to  be 
complex.  The  source  and  field  points  are  not  restricted  to  the 
water  column;  one  or  the  other,  or  both,  may  be  within  the 
seafloor.  Expression  (1)  assumes  that  the  seafloor  behaves  as 
an  acoustic  fluid,  thus  the  model  does  not  include  shear 
waves.  In  the  remainder  of  this  article,  the  scattered  pressure 
field  at  position  r,  due  to  a  unit  point  source  (defined  in  the 
preceding  text)  at  position  r0  will  be  denoted  ps( r,  r0)  and 
defined  as  the  first-order  solution  for  the  field.  To  simplify 
the  exposition,  essential  results  will  be  presented  in  this  sec¬ 
tion  with  the  derivation  deferred  to  Appendix  A.  The  deriva¬ 
tion  follows  a  method  given  by  Ivakin28  and  requires 
somewhat  abstract  notation  for  the  general  result.  To  avoid 
immediate  involvement  in  notational  issues,  the  discussion 
will  be  initiated  with  a  less  general  problem:  semi-infinite 
water  and  sediment  half-spaces,  each  being  homogeneous 
with  regard  to  sound  speed  and  density.  The  boundary 
between  these  two  media  will  be  assumed  to  be  flat  except 
for  the  small-scale  roughness.  The  scattered  field  at  a  point 
r  =  (x,  y,  z)  in  the  water  column  due  to  a  unit  source  in  the 
water  column  at  r0  =  (x0,  yo>  zo )  is  given  by  the  expression 


Mr,ra)  =  ^ 


co2(K2-Kl)G(0)(r,,r)G(0)(r,,r0) 

V,G^(r,,r)|z'=0+ 

(2) 


dx  dy 

1  1 
Pi  Pi 
V/G(°)(r',ro)|z,=0-lc(V,y) 


The  integration  is  over  the  horizontal  plane  at  z!  =  0  that  sep¬ 
arates  water  and  seafloor,  thus,  r'  =  (x',  /,  0).  The  perturbed 
water- seafloor  boundary  is  given  by  z!  =  £(x',  /).  The  sub¬ 
scripts  1  and  2  refer  to  water  and  seafloor,  respectively,  and 
the  corresponding  compressibilities  are  denoted 

Ki  =  ,  i  =  1,2  (3) 

CJPi 

with  ct  being  the  sound  speed  and  pt  being  the  density. 
Attenuation  in  the  seafloor  is  incorporated  by  allowing  k2 
(hence  c2)  to  be  complex. 

In  Eq.  (2),  the  Green’s  functions  for  the  zeroth-order 
problem  (two  semi-infinite  half  spaces  separated  by  the  plane 
z!  =  0)  appear.  The  zeroth-order  solution  for  this  case  can  be 
found,  e.g.,  by  means  of  wavenumber  integration.  The  two 
Green’s  functions  of  interest  are  G(0)(r',  r0)  for  the  field  at 
the  integration  point  r'  due  to  a  source  at  the  true  source  posi¬ 
tion  r0,  and  G(0)(r',  r)  for  the  field  at  the  integration  point  due 
to  a  source  at  the  field  point  r  The  gradient  of  G(0)(r',  r0)  is  to 
be  evaluated  immediately  below  the  interface,  while  the 
gradient  of  G(0)(r',  r)  is  to  be  evaluated  immediately  above 
the  interface.  As  the  pressure  field  is  continuous  across  the 
interface,  the  factor  G(0)(r',  r)G(0)(r',  r0)  can  be  evaluated 
either  above  or  below.  In  fact,  the  same  is  true  for  the  x'-  and 
y' -partial  derivatives.  The  normal  gradient  of  the  pressure 
field,  however,  is  discontinuous  across  the  interface,  so  care 
must  be  taken  in  the  evaluation  of  the  z!  component  of  the 
gradient.  As  noted  in  Appendix  A,  the  displacement  matching 
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condition  that  relates  the  normal  derivative  values  on  either 
side  of  the  interface  can  be  used  to  show  that  it  does  not  mat¬ 
ter  which  of  the  two  normal  derivatives  in  the  dot  product  of 
Eq.  (2)  is  evaluated  above  as  long  as  the  other  is  evaluated 
below.  Using  this  symmetry,  it  follows  that  the  scattered  field 
obeys  reciprocity,  that  is,  ps(r ,  r0)=p5(r0,  r).  An  equally 
important  property  of  this  solution  is  that  it  agrees  with  con¬ 
ventional  small-roughness  perturbation  theory.  In  particular, 
if  the  mean-square  scattered  pressure  is  used  to  extract  a 
scattering  cross  section,  the  result  will  agree  with  formally- 
average  perturbation  theory  as  will  be  shown  in  the  next 
section. 

Expression  (2)  applies  when  both  the  source  and  field 
points  are  in  the  water  column,  z0  >  0,  z  >  0.  The  case  z0  >  0, 
z  <  0  is  of  interest  in  modeling  the  acoustic  field  that  pene¬ 
trates  the  seafloor.  This  case  has  been  published  by  Thorsos 
and  co-workers29  with  a  resulting  expression  that  is  equiva¬ 
lent  to  Eq.  (2)  except  with  the  density  prefactor  px  replaced 
by  p2 •  The  opposite  case,  with  source  in  the  seafloor  and 
field  point  in  the  water  gives  an  expression  identical  to 
Eq.  (2),  including  the  prefactor  px.  The  remaining  case,  with 
both  source  and  field  point  in  the  seafloor  (zo<0,  z<0), 
uses  the  same  expression  as  that  for  source  in  the  water  and 
field  point  in  the  sediment.  Note  that  identity  of  expressions 
does  not  imply  identity  of  scattered  pressures  because  the 
unperturbed  Green’s  functions  that  enter  these  expressions 
may  have  different  values.  All  of  the  preceding  relationships 
lead  to  the  reciprocity  statement  ps(rb,  ra)/pb  =ps(ra,  r b)/pa 
given  in  Ref.  21  (Chap.  8)  for  the  zeroth-order  case,  yet  ap¬ 
plicable  in  all  cases. 

Expression  (2)  is,  in  fact,  much  more  general  than  might 
be  first  supposed.  It  is  applicable  to  horizontally  stratified 
waveguides  if  the  zeroth-order  Green’s  functions  are  eval¬ 
uated  for  the  waveguide.  The  source  and  field  points  may  be 
either  in  the  water  or  seafloor,  and  the  density  prefactor  must 
be  evaluated  at  the  field  point.  The  densities,  pt,  and  com¬ 
pressibilities,  Kh  in  the  integrand  are  to  be  evaluated  at 
z'  =  0+  if  i  =  1  (water)  or  z!  =  0“  if  /  =  2  (seafloor).  Thus  the 
present  model  can  be  applied  to  range-independent  wave¬ 
guides,  where  modal  calculations  can  be  used  to  evaluate  the 
zeroth-order  Green’s  functions  and  their  gradients.  The  fol¬ 
lowing  section  contains  an  example  in  which  this  model  is 
applied  to  a  finely  layered  seafloor. 

All  of  the  cases  discussed  in  the  preceding  text  are  spe¬ 
cial  cases  of  a  more  general  result  for  range-dependent 
waveguides,  derived  in  Appendix  A.  The  general  setting  is 
depicted  in  Fig.  1,  and  the  scattered  field  is  given  by 


P*(r,r0)  =  [®2(k2  -  Ki)]G(0)(r',r)G(0)(r',r0) 


_(l_l)V'G(0)(r',r)|lv,=o+ 
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•V'G(V,ro)  Uo-1  C(r'). 
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Because  the  zeroth-order  water- seafloor  interface  is  not  flat, 
a  surface  integral  replaces  the  .^-/-integral  in  Eq.  (2).  The 
coordinate  w'  is  normal  to  the  zeroth-order  interface  at  every 
point.  Symmetry  with  respect  to  evaluation  of  the  two 


Green’s  function  factors  at  either  w' =  0+  or  wf  =  0  still 
applies,  so  reciprocity  is  respected.  The  densities,  ph  and 
compressibilities,  Kh  are  to  be  evaluated  at  w'  =  0+  if  /  =  1 
(water)  or  w'  =  0“  if  /  =  2  (seafloor).  Note  that  these  parame¬ 
ters  may  be  functions  of  the  transverse  coordinates  involved 
in  the  surface  integration  and  that  the  compressibility  is  a 
complex  function  of  position  in  the  seafloor.  Furthermore,  the 
scattering  interface  is  not  restricted  to  the  water-seafloor  inter¬ 
face;  it  can  be  a  buried  interface  in  the  sediment,  a  sediment- 
basement  interface,  for  example.  Finally,  the  small  roughness 
£(r)  is  always  measured  normal  to  the  local  tangent  plane  of 
the  slowly  varying  surface. 

Computation  of  the  Green’s  functions  and  their  gradients 
is  more  difficult  in  the  range-dependent  case.  In  this  article,  a 
parabolic  equation  method  will  be  used  with  special  effort  to 
compute  the  gradients  near  the  zeroth-order  interface.  In 
Sec.  V,  reverberation  in  a  wedge-shaped  waveguide  will  be 
computed  as  an  example  of  the  application  of  Eq.  (4). 

III.  BOTTOM  BACKSCATTER  FOR  HOMOGENEOUS 
HALF-SPACE 

Equation  (2)  of  the  previous  section  will  be  applied  to  a 
simple  backscatter  scenario  for  which  the  theoretical  scatter¬ 
ing  cross  section  is  available.  This  provides  a  partial  valida¬ 
tion  of  the  simulation  results.  The  environment  is  appropriate 
for  bottom  backscatter  measurement  where  the  sea  surface  is 
far  away  so  its  effects  can  be  easily  separated  in  time.  As 
such,  the  water  depth  is  assumed  to  be  infinite,  and  the  me¬ 
dium  consists  of  two  fluid  half-spaces  separated  by  a  rough 
seafloor.  The  upper  half- space  is  homogeneous  water  with 
sound  speed  cx  and  density  px.  The  attenuation  in  water  is 
ignored.  The  lower  half-space  is  homogeneous  sediment  with 
sound  phase  speed  c2p ,  density  p2 ,  and  attenuation  factor  S2. 
The  complex  speed  is  c2  =  c2p/(l  +  id2).  The  coordinate  sys¬ 
tem  is  defined  such  that  the  x  and  y  axes  are  horizontal.  The 
random  seafloor  is  given  by  z  =  £(x,  y),  assumed  to  have  zero 
mean.  The  z  axis  points  upward  with  the  mean  seafloor  as  its 
origin.  This  is  a  special  case  of  the  general  rough  interface 
discussed  in  the  previous  section  because  the  slowly  varying 
component  of  the  roughness  is  simply  a  constant. 

The  source  and  receiver  are  co-located  at  (0,0,zo).  The 
backscattered  sound  field  is  to  be  calculated  for  a  number  of 
independent  realizations  of  the  rough  seafloor.  These  calcu¬ 
lations  will  result  in  a  set  of  backscattered  time  series,  which 
will  form  the  statistical  ensemble  for  estimating  the  bottom 
scattering  cross  section.  The  theoretical  scattering  cross  sec¬ 
tion  for  this  case  is  known21  (Chap.  13)  and  will  be  com¬ 
pared  with  our  simulation  results.  This  section  contains 
many  details  that  will  be  common  for  other  sections.  It  is  di¬ 
vided  in  subsections  to  make  it  easier  to  follow. 

A.  Synthesis  of  rough  surfaces 

In  this  special  backscatter  geometry  with  range- 
independent  environment,  the  only  quantity  depending  on  the 
azimuth  angle  (j)  under  the  integral  (2)  is  the  rough  surface 
£(x,  y)  =  £(R,  (j))  when  the  horizontal  dimensions  are  ex¬ 
pressed  in  polar  coordinates.  Therefore,  if  one-dimensional 
realizations  of  the  roughness  are  obtained  by  pre-summing 
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over  the  azimuth,  rj(R)  =  j^n  £(/?,  </>)d(/>,  Eq.  (2)  is  simplified 
to  the  one-dimensional  integral: 


W(  K)  = 
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where  K  =[KX,  Ky\ ,  K0  =  2n/L0 ,  and  L0  defines  the  outer 
scale  of  the  roughness.  For  examples  in  this  section  and  next, 
y2  =  3.25,  L0  =  10  m,  and  w2  =  8.496  x  10-6  m4-y2.  The  pro¬ 
cedure  of  obtaining  realizations  of  rj(R)  =  £(R,  (j))d(j) 

through  Hankel  transforms  is  given  in  Appendix  B  and  can 
also  be  found  in  Li  et  al .30 

B.  Zeroth-order  field 


where  the  matching  relationship  between  the  vertical  gra¬ 
dients  in  the  upper  and  lower  media  has  been  used.  For  con¬ 
venience,  the  dependence  of  the  Green’s  functions  on  the 
source  and  receiver  coordinates  [both  equal  to  (0,0, z0)]  has 
been  suppressed.  The  characteristics  of  the  rough  interface 
C(x,  y )  are  specified  by  the  two-dimensional  power-law 
power  spectrum 


To  calculate  the  backscattered  held  in  Eq.  (5),  the  zeroth- 
order  pressure  held  and  its  horizontal  and  vertical  derivatives 
on  the  hat  interface  are  needed.  For  the  present  range- 
independent  problem,  these  quantities  can  be  obtained  by  sep¬ 
aration  of  variables  and  wavenumber  synthesis31  (Chap.  1). 
They  are 
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where  V  is  the  plane-wave  rehection  coefficient: 


P2P1  ~  Pl/^2 

P2P1  +  Pi  P2 


(8) 


with  Pi  =  yjk\—  K2,  and  /?2  =  \/k2—  K2.  In  the  hrst  equa¬ 
tion  of  (7),  the  hrst  term  is  the  direct  arrival  from  the  source 
to  the  seafloor  and  the  second  term  the  reflected  arrival.  The 
quantities  in  Eq.  (7)  are  numerically  evaluated  along  the 
interface  to  be  fed  into  Eq.  (5).  When  the  integration  point  is 
in  the  far  held  of  the  source/receiver  (kiz0^>  1),  a  ray  form 
of  Eq.  (7)  can  be  used  through  the  stationary  phase 
approximation: 
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where  =  k\  sin  G),  //{ =  \J k\  ~  cosG))2,  tan  G0  =  zo /R' , 
and  V  0  is  the  rehection  coefficient  evaluated  at  the  grazing 
angle  0O-  When  these  approximate  zeroth-order  terms  are 
used,  the  scattered  held  is  simplihed  to 
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where  ak  =  k2/ki.  It  is  shown  in  Appendix  D  that  Eq.  (10) 
results  in  the  same  expression  for  backscattering  cross  sec¬ 
tion  as  predicted  by  hrst-order  perturbation  theory.  When  the 
zeroth-order  helds  are  obtained  through  either  Eq.  (7)  or  (9), 
the  hrst-order  scattered  held  at  any  frequency  can  be  calcu¬ 
lated  using  Eq.  (5)  or  (10).  For  backscatter  simulations,  the 
zeroth-order  held  from  the  source  to  the  scatterers  and  that 
from  the  scatterers  to  the  receiver  are  the  same,  hence  the 
computation  time  is  half  of  that  for  the  bistatic  geometry. 
For  different  realizations  for  the  rough  interface,  the  same 
zeroth-order  helds  are  used,  making  it  efficient  to  calculate 
large  numbers  of  scattered  held  realizations. 


C.  Fourier  synthesis 

Fourier  synthesis  is  used  to  obtain  a  simulated  pressure 
time  series.  In  this  case,  the  elapsed  time  after  transmission 
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provides  angular  resolution  of  the  scattered  signal,  hence  of 
the  scattering  cross  section.  A  Gaussian  pulse  with  center 
frequency  fc  is  chosen  in  the  form 

s(t)  =  se(t)  cos(27 lfct),  (12) 

where 

se(t)  =  exp[— 2(7iQ/)2]  (13) 

is  the  baseband  envelope.  The  pulse  peaks  at  time  zero  with 
value  unity  and  falls  off  exponentially  versus  the  square  of 
time.  The  baseband  envelope  has  the  following  Fourier 
transform: 

Seif)  =  T  Jse(f)  exp  {ilnft))dt 

=  — |p=exp(— (14) 

Q\/2tt  V  2Q7 

The  bandwidth  is  controlled  by  the  parameter 

Q  =fb/v/iogJ?),  (15) 

where  fb  is  the  baseband  frequency  at  which  the  signal  fre¬ 
quency  spectrum  equals  half  peak  power  (fb  =fc/ 40  in  our 
simulations,  giving  a  bandwidth  of  5%  of  the  center  fre¬ 
quency).  The  highest  difference  frequency  on  either  side  of 
the  center  frequency  is  chosen  to  be  at  a  value  where  the 
spectrum  is  X  decibels  below  the  peak  spectrum, 

=  <16> 

The  value  X  =  60  is  used  for  all  examples  in  this  paper.  This 
determines  the  bandwidth  over  which  the  scattered  field 
must  be  calculated.  The  choice  of  the  frequency  increment  is 
dictated  by  the  time  window  over  which  the  time-domain 
scattered  field  is  desired.  This  is  achieved  by  choosing  a 
maximum  range  rmax  at  which  the  roughness  terminates. 
This  termination  point  determines  the  smallest  grazing  angle 
that  can  be  obtained. 

D.  Results 

For  the  set  of  parameters  given  in  Tables  I  and  II, 
the  zeroth-and  first-order  scattered  field  for  a  set  of  100 


TABLE  I.  Environmental  parameters  used  in  simulations.  All  simulations 
use  water  sound  speed  1500  m/s  and  density  1000  kg  m-3. 


Case 

Layer  Thickness  Sound  speed 

Density 

Attenuation 

No. 

(m) 

(m/s) 

(kg  m-3) 

(dB/A) 

Two  half  spaces 

1 

oo 

1700 

2000 

0.5 

Approximate  Lambert 

1 

0.2 

1446.2 

1400 

0.05 

2 

oo 

1700 

2000 

0.5 

Waveguide 

1 

oo 

1700 

2000 

0.5 

Wedge 

1 

oo 

1700 

1500 

0.5 

realizations  of  the  rough  seafloor  is  obtained.  The  roughness 
relief  realizations  are  statistically  independent.  Conse¬ 
quently,  the  scattered  intensity  realizations  are  also  statisti¬ 
cally  independent.  Thus  the  variance  of  the  mean  intensity 
from  N  realizations  is  reduced  by  a  factor  of  1  /N  compared 
to  a  single  realization.  For  example,  100  realizations  will 
reduce  the  variance  to  1/100,  or  the  standard  deviation  to 
1/10  of  the  mean. 

Small-roughness  perturbation  theory  should  be  accurate 
for  this  example  in  which  the  product  of  RMS  roughness  and 
acoustic  wavenumber  is  0.115.  The  numerical  data  are  used 
to  estimate  the  backscattering  cross  section  following  the 
steps  given  in  Appendix  C.  The  results  are  compared  to  the 
theoretical  backscattering  cross  section  in  Fig.  2.  The  model 
is  for  the  same  environment  as  the  simulations,  and  the  fre¬ 
quency  is  1000  Hz.  The  angular  range  shown  for  the  compar¬ 
ison  is  4°  to  60°,  over  which  the  two  are  consistent, 
including  the  abrupt  change  near  28°,  the  critical  angle. 

IV.  APPROXIMATE  LAMBERTIAN  SCATTERING 

Lambert’s  law  is  a  bottom  scattering  cross  section 
widely  used  in  both  applied  and  theoretical  models  because 
of  its  simplicity.  For  backscatter,  the  cross  section  is  propor¬ 
tional  to  the  second  power  of  the  sine  of  the  grazing  angle: 

cr/(0)  =  jusin2  0,  (17) 

where  /i  is  a  constant.  Lambert’s  law  is  a  concept  related  to 
multiple  scattering  in  low-loss  media,  and  single- scattering 
theory  cannot  give  Lambert’s  law  over  all  angles  and  fre¬ 
quency.  However,  it  is  desirable  to  have  a  physics-based 
model  that  resembles  the  Lambert  model  in  some  useful 
angular  and  frequency  range,  so  that  results  from  the 
physics-based  model  can  be  compared  to  those  obtained 
from  codes  that  assume  Lambertian  scattering.  Here  first- 
order  rough  interface  scattering  theory  is  applied  to  mimic 
Lambertian  scattering.  Because  small  grazing  angles  are  of¬ 
ten  the  important  regime  for  long-range  propagation  and 
reverberation,  our  limited  goal  is  to  find  a  region  of  small 
grazing  angles  for  a  given  frequency  where  the  first-order 
rough  interface  theory  fits  Lambert’s  law. 

As  the  grazing  angle  approaches  zero,  the  interface  scat¬ 
tering  cross  section  also  necessarily  goes  to  zero  as  required 
by  conservation  of  energy.  The  scattering  cross  section  of 
first-order  rough  surface  scattering  theory  approaches  zero  as 
the  fourth  power  of  grazing  angle,  the  same  as  that  predicted 
for  the  rough  sea  surface  (pressure  release  boundary  condi¬ 
tion).  For  this  theory  to  mimic  Lambertian  scattering,  which 
is  proportional  to  the  second  power  of  the  grazing  angle,  a 
thin,  soft  layer  above  the  rough  interface  is  added,  hence 
making  the  local  grazing  angle  steeper.  Adding  a  single  layer 
can  only  fit  the  Lambert  cross  section  in  a  finite  range  of 
grazing  angles  for  fixed  frequency.  Adding  multiple  layers, 
in  principle,  can  increase  the  angular  range  of  fit.  In  adding  a 
single  layer,  one  new  feature  of  the  theory  will  also  be 
shown. 

The  environment  is  shown  in  Fig.  3,  where  all  the  envi¬ 
ronmental  parameters  are  the  same  as  in  those  of  the 
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TABLE  II.  Geometrical  and  other  simulation  parameters. 


Case 

fc  (Hz) 

Source  coordinate  (m) 

Rec.  coordinate  (m) 

No.  Realiz. 

Surface  length,  L  (m) 

No.  points,  N 

Two  half  spaces 

1000 

150  (From  bottom) 

150  (From  bottom) 

100 

8,000 

20,000 

Approximate  Lambert 

1000 

150  (From  bottom) 

150  (From  bottom) 

100 

8,000 

20,000 

Waveguide 

250 

15  (From  surface 

25  (From  surface) 

100 

10,000 

10,000 

Wedge 

250 

35  (From  surface) 

35  (From  surface) 

100 

4,000 

4,000 

previous  section  except  that  a  20-cm  mud  layer  with  sound 
speed  lower  than  the  water  sound  speed  is  put  above  the 
rough  sand  sediment  half-space  (Table  I).  The  thin,  slow 
layer  will  tend  to  increase  the  local  grazing  angle  when  the 
incident  acoustic  wave  impinges  on  the  rough  interface, 
hence  decreasing  the  rate  at  which  the  cross  section 
approaches  zero  at  small  grazing  angles. 

The  formula  applicable  to  this  problems  is 


Ps(z o)  = 


P(z  o) 

471 


»oc 

R'dR'r]{R') 

o 


(02(k2  -  Ki)(G(0)(r'))2 


j_\  /dG°\2 
Pi)  \dR'J 


Plfl 

pApi 


/0G<°V)\ 
V  dzf  ) 


(18) 


This  expression  is  very  similar  to  Eq.  (5)  with  the  following 
differences:  First,  the  mean  of  the  rough  interface  is  at 
z  =  —0.2  m  instead  of  zero.  Second,  the  density  pi  =  1400 
kg/m3  refers  to  the  density  in  the  mud  layer,  rather  than  that 
in  water.  The  density  p2  remains  that  for  the  sand  half- space 
(Table  I).  Third,  the  Green’s  function  is  one  linking  the 
source  to  a  point  on  the  mean  interface  at  z  =  —0.2  m.  The 
Green’s  function  and  its  gradients  are  numerically  calculated 
using  wavenumber  integration  for  the  three-layer  scenario. 
All  other  quantities  are  the  same  as  in  the  previous  section, 
and  simulation  parameters  are  given  in  Table  II.  The  results 
in  the  form  of  backscattering  strength  are  shown  in  Fig.  4, 


FIG.  2.  Comparison  of  backscattering  strength  as  a  function  of  grazing 
angle  calculated  from  the  mean-square  pressure  of  100  realizations  to  the 
perturbation  theory  model  for  the  two-half-space  case. 


where  the  simulation  agrees  with  perturbation  theory  as 
expected.  Both,  in  turn,  agree  with  Fambert’s  law  over  the 
angle  range  of  5°  to  28°  with  the  latter  being  the  critical 
angle.  This  demonstrates  that  one  can  use  the  three-layer 
model  to  mimic  Fambertian  scattering  over  these  angles  in 
the  neighborhood  of  1  kHz. 


V.  REVERBERATION  IN  2D  WAVEGUIDES 

The  first-order  perturbation  result  will  be  applied  to  two 
waveguide  reverberation  problems.  The  first  is  a  range- 
independent  problem  for  which  a  solution  from  numerically 
averaged  perturbation  theory  is  available  for  comparison. 
The  second  is  reverberation  in  a  wedge-shaped  waveguide. 
Both  examples  are  in  two-dimensional  geometry  where  the 
environment  has  no  y  dependence.  In  both  cases,  the  para¬ 
bolic  equation32  is  used  to  calculate  the  zeroth-order  Green’s 
function.  The  use  of  PE  to  calculate  reverberation  has  been 
reported  previously  by  Fingevitch  and  FePage.27  The  differ¬ 
ence  here  is  that  the  parabolic  equation  is  used  to  calculate 
the  Green’s  function  corresponding  to  the  slowly  changing 
environment,  leaving  the  scattering  problem  to  be  handled 
by  perturbation  theory,  hence  saving  computation  time. 

The  first  problem  is  one  defined  by  the  Reverberation 
Modeling  Workshops  sponsored  by  the  U.S.  Navy.9,33  The 
environment  is  a  50-m  deep  isovelocity  waveguide  with  a  ho¬ 
mogeneous  sandy  sediment  half-space.  The  water- sediment 
interface  is  rough  with  the  following  ID  power  spectrum  that 
was  defined  for  the  “typical”  roughness  case:9,33 

w'°iK)  =  Afhoy  <19) 

where  /z  =  0.3 16  m,  and  KL  =  0.1  m-1  for  the  numerical 
examples  for  the  rest  of  the  paper.  The  source  depth  is  15  m 


Water 

c  =  1 500  m/s 

\ 

speed  ratio  =  0.964,  density  ratio  =  1 .4,  atten  =  0.05  dB/wavelength  0.2  m 

speed  ratio  =  1.133,  density  ratio  =  2.0,  atten  =  0.5  dB/wavelength 


FIG.  3.  Layered  seafloor  model  used  to  mimic  Lambertian  scattering  over 
limited  angular  and  frequency  ranges. 
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6  (deg) 


FIG.  5.  Comparison  of  mean  reverberation  level  vs  time  from  100  rough 
surface  realizations  using  the  PE  and  normal  mode  methods. 


FIG.  4.  Comparison  of  model  and  simulation  approximating  Lambert’s  law. 
The  simulation  curve  and  the  model  curve  from  perturbation  theory  are 
consistent  for  all  grazing  angles.  These  two  are  consistent  with  Lambert’s 
law  over  the  angle  range  of  5°  to  28°. 


and  the  receiver  depth  is  25  m  at  the  same  horizontal  range. 
The  narrow-band  source  has  a  center  frequency  of  250  Hz. 
Details  of  all  parameters  involved  are  found  in  Yang  et  al9 
The  formula  for  this  problem  takes  on  the  following 
form: 


|*oc  [—  J  r 

ps(z,  z0)=-1  V)  ®2(k2-ki) 

471 J  o 

x  (g(°\z0,x' ,z')G(°\z,x' 


1  1 

P2  P\ 

/9G^(zo,t/,z/)9G°(z,x/,z/)\  p2  (  1  1 


\  dx'  dx' 

fdG^(z0,x',z')  dG<®(z,jd,z!) 


Pi  \P 2  Pi 


V  dzf 


dz' 


(20) 


In  this  expression,  z0  is  the  source  depth,  and  z  is  the  receiver 
depth.  The  Green’s  functions  connect  the  source  at  (0,  z0)  or 
the  receiver  at  (0,  z)  to  points  just  above  the  mean  surface 
(x\  z').  These  Green’s  functions  and  their  gradients  for  the 
range-independent  problem  can  be  obtained  through  a  num¬ 
ber  of  methods,  including  the  normal  mode  method.  The  PE 
method  was  chosen  here  because  the  PE  approach  can  also 
be  used  for  range-dependent  problems.  Once  the  Green’s 
functions  and  their  gradients  are  obtained,  the  same  steps  as 
those  in  Sec.  Ill  can  be  taken  to  calculate  reverberation  time 
series  for  different  rough  surface  realizations.  Simulation  pa¬ 
rameters  are  given  in  Table  II.  Figure  5  compares  the  results 
using  the  PE  method  to  those  obtained  using  the  normal 
mode  method.9  The  two  results  are  essentially  identical  for 
times  greater  than  1  s.  The  discrepancy  between  the  two  for 
times  shorter  than  1  s  is  due  to  the  fact  that  neither  method 
treats  high-wavenumber  propagation  properly  as  both  are 
intended  only  for  computation  of  long-range  reverberation. 
The  purpose  of  the  comparison  is  primarily  to  establish 
using  the  PE  as  a  viable  method  for  range-dependent  rever¬ 
beration  problems.  Incidentally,  favorable  results  from  the 


energy-flux18  and  ray-tracing  methods  for  this  range- 
independent  example  can  be  found  in  Ref.  33. 

Reverberation  in  a  wedge-shaped  waveguide  is  dis¬ 
cussed  next.  The  waveguide  is  the  “ASA  wedge”  (Ref.  34). 
The  waveguide  has  iso  velocity  sound  speed  of  1500  m/s  and 
a  sloping  bottom  of  2.86°.  The  bottom  sound  speed  is  1700 
m/s  with  an  attenuation  coefficient  of  0.5  dB  per  wavelength 
and  sediment-to-water  density  ratio  of  1.5  (Table  I).  The 
source  and  receiver  are  co-located  at  35  m  depth  where 
the  water  depth  is  200  m.  The  acoustic  source  spectrum  is 
the  same  as  in  the  range-independent  waveguide.  The  for¬ 
mula  for  this  case  is 


Kl)(G(°)(zo/))2 


2_V— V 

Pi)  V  dl  ) 


inn 

Pl\P2 


dG^v 

dn'  ) 


above  bottom 


(21) 


Several  items  in  this  expression  need  clarification.  First,  /'  is 
along  the  sloping  bottom.  Second,  the  rough  interface  £(/)  is 
measured  normal  to  the  sloping  bottom.  Third,  the  gradient 
in  the  last  term  involving  dn'  is  the  gradient  normal  to  the 
sloping  bottom.  Finally,  the  Green’s  function  and  its  gra¬ 
dients  are  calculated  using  PE  and  evaluated  along  and  just 
above  the  sloping  bottom.  Because  the  source  and  receiver 
are  co-located,  the  Green’s  function  needs  only  to  be  calcu¬ 
lated  once  for  each  frequency,  and  the  quantities  needed  are 
the  pressure  field  and  its  horizontal  and  normal  gradients 
along  the  water-bottom  interface.  If  the  source  and  receiver 
are  not  co-located,  then  the  Green’s  functions  from  both  the 
source  and  the  receiver  to  the  water-bottom  interface  will  be 
required,  making  the  computational  time  twice  as  long  as 
that  for  the  monostatic  case.  The  number  of  frequencies  for 
which  the  Green’s  functions  are  needed  depends  on  the 
acoustic  source  bandwidth.  The  number  of  realizations  does 
not  impact  the  computational  load  on  the  Green’s  functions. 

The  reverberation  level  is  calculated  from  100  realiza¬ 
tions  with  the  simulation  parameters  given  in  Table  II  and 
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FIG.  6.  Mean  reverberation  level  vs  time  from  100  realizations  for  a  wedge- 
shaped  waveguide. 


with  the  results  given  in  Fig.  6.  One  interesting  feature  is  that 
the  reverberation  level  remains  roughly  a  constant  for  time 
greater  than  1  s.  It  is  known  that  sound  propagating  upslope 
tends  to  gradually  shift  energy  to  higher  grazing  angles.  The 
portion  of  energy  shifted  to  angles  above  the  critical  angle 
will  be  lost  in  the  sediment.  From  a  modal  point  of  view, 
high-order  modes  will  be  cut  off  first,  followed  by  lower- 
order  modes.  Therefore,  there  seems  to  be  a  contradiction 
that  while  the  total  propagating  energy  decreases  going 
upslope,  the  reverberation  level  does  not  decrease.  A  closer 
examination  reveals  that  the  sound  intensity  along  the  sloping 
bottom  does  not  have  a  strong  decrease  until  all  modes  are 
cut  off  toward  the  end  of  the  wedge.  Because  the  scatterers 
are  on  the  bottom,  the  nearly  range-independent  intensity 
along  the  bottom  is  consistent  with  the  observed  steady  rever¬ 
beration  level  for  the  two-dimensional  case. 

VI.  SUMMARY 

A  rough-interface  reverberation  model  has  been  devel¬ 
oped  using  first-order  perturbation  theory,  where  the  unper¬ 
turbed  background  interface  is  allowed  to  be  arbitrary  in 
shape,  and  the  medium  surrounding  the  rough  interface  can 
be  heterogeneous  and  layered.  The  advantage  of  this 
approach  for  calculating  the  reverberation  field  is  that  the 
slowly  changing  unperturbed  environment,  which  does  not 
contribute  to  backscatter  but  alters  forward  propagation,  can 
be  handled  by  fast  numerical  methods  such  as  the  parabolic 
equation.  In  this  model,  the  reverberation  field  is  calculated 
for  realizations  of  the  small  roughness  superimposed  on  the 
slowly  varying  interface.  The  overall  increased  computa¬ 
tional  speed  comes  from  the  fact  that  numerical  steps  in 
range  need  only  accommodate  the  slow  range  dependence 
and  are  much  greater  than  the  acoustic  wavelength.  This 
contrasts  with  the  wavelength-scale  range  steps  required  by 
Lingevitch  and  LePage,27  for  example.  The  computation 
speed  is  further  enhanced  because  the  forward  propagating 
field  is  common  to  all  realizations;  obtaining  large  numbers 
of  realizations  does  not  require  additional  time  proportional 
to  the  number  of  realizations.  The  primary  purpose  of  this 
model  is  to  provide  a  physics-based  approach  that  is  reason¬ 
ably  fast,  so  the  accuracy  of  faster,  but  possibly  less  accu¬ 
rate,  models  can  be  checked. 
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APPENDIX  A:  PERTURBATION  THEORY 

This  appendix  provides  a  derivation  of  the  main 
perturbation-theoretic  result,  (4).  The  derivation  employs  the 
unified  approach  given  by  Ivakin28  in  which  roughness  scat¬ 
tering  is  treated  as  a  form  of  volume  scattering.  Deformation 
of  an  interface  causes  a  position-dependent  change  in  acous¬ 
tic  parameters,  and  an  exact  integral  equation  is  developed 
for  the  resultant  change  in  the  pressure  field.  This  equation 
can  be  expressed  in  the  following  form: 

G(r,r0)  =  G(0)(r,r0) 

+  ^|  {[®2(k  -  Ko)]G(0)(r',  r)G(r',  r0) 

I-3-W'G(<V,r)  •  V'G( r',ro)}<fV.  (Al) 
P  Po)  J 

As  in  all  other  parts  of  this  article,  exp(-/cof)  time  depend¬ 
ence  is  assumed.  Equation  (Al)  includes  a  prefactor  p(r)  not 
appearing  in  Ref.  28  owing  to  differing  definitions  of  the 
Green’s  function.  This  expression  is  quite  general  with  no 
limitations  on  the  structure  of  the  acoustic  media  (apart  from 
the  assumption  of  fluid  behavior)  and  with  no  restriction  on 
the  location  of  the  source  and  receiver.  With  reference  to 
Fig.  1,  G(0)(r,  r0)  in  Eq.  (Al)  is  the  “zeroth-order”  Green’s 
function  for  the  unperturbed,  smooth  bathymetry.  It  gives 
the  field  at  point  r  due  to  a  unit  point  source  situated  at  r0. 
This  Green’s  function  is  the  solution  of 


Po(r)V 


•  — j-— VG^(r,  r0) 

LPoM 

4n5(r  -  r0), 


+  kl(r)G{0)(r,  r0) 


(A2) 


where  p0(r)  is  the  position-dependent  mass  density  (for  both 
the  unperturbed  seafloor  and  water  column),  and  £0(r)  is  the 
position-dependent  unperturbed  wavenumber.  Wavenumber 
and  sound  speed  c  are  related  by  co  =  ck ,  and  both  sound 
speed  and  wavenumber  may  be  complex  with  imaginary 
parts  that  determine  attenuation.  In  Eq.  (Al),  k  is  the 
position-dependent  compressibility,  related  to  density  and 
sound  speed  as  follows: 

k=4~-  (a3) 

czp 

The  position  dependencies  of  compressibility,  sound  speed, 
and  density  have  not  been  indicated  in  Eq.  (A3)  for  the  sake 
of  simplicity.  This  equation  refers  to  the  true  seafloor  acous¬ 
tic  parameters,  but  the  same  relation  holds  for  the  unper¬ 
turbed  parameters. 

The  “true”  Green’s  function  in  Eq.  (Al)  is  denoted 
G(r,  r0),  and  is  the  solution  of  Eq.  (1)  in  which  the  density 
and  wavenumber  are  given  their  true  values,  that  is,  the  val¬ 
ues  corresponding  to  the  true  bathymetry  indicated  in  Fig.  1. 
Departure  of  the  true  bathymetry  from  the  smooth 
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bathymetry  (solid  line  in  Fig.  1)  results  in  changes  in  density 
and  compressibility  with  the  new  position-dependent  acous¬ 
tic  parameters  indicated  in  Eq.  (Al)  by  the  symbols  p  and  k. 
It  is  important  to  realize  that  the  new  parameters  differ  from 
the  old  only  in  small  lens- shaped  regions  near  the  smooth 
interface.  In  particular,  for  the  water-sediment  interface,  the 
parameters  change  from  water  to  sediment  values  for  those 
regions  where  the  true  interface  lies  above  the  smooth  inter¬ 
face.  Similarly,  they  change  from  sediment  to  water  values 
for  those  regions  where  the  true  interface  lies  below  the 
smooth  interface.  For  buried  interfaces,  the  changes  are 
those  dictated  by  the  two  different  sediment  types. 

Use  of  the  term  “perturbed”  suggests  that  the  changes 
between  the  parameters  for  the  true  and  zeroth-order  prob¬ 
lems  are  small.  In  fact,  Eq.  (Al)  is  exact,  and  these 


perturbations  can  be  arbitrarily  large.  In  this  article,  how¬ 
ever,  the  change  in  bathymetry  is  assumed  to  be  small,  so 
that  a  relatively  simple  perturbative  solution  of  Eq.  (Al)  is 
possible.  This  solution  method  parallels  one  given  by 
Ivakin.28 

The  scattered  field  will  be  defined  as  the  difference 
between  the  zeroth-order  field  and  the  exact  field.  In  the 
small-roughness  approximation,  the  scattered  field  is 
computed  to  first  order  in  the  perturbation  of  the  bathym¬ 
etry,  here  denoted  £(r'),  and  illustrated  in  Fig.  7.  The 
volume  integration  in  Eq.  (Al)  is  taken  over  the  thin 
lens- shaped  regions  between  the  perturbed  and  unper¬ 
turbed  interfaces.  The  scattered  field  to  first  order  in  £(r') 
can  be  written  as 


Psix,  r0)  =  ^fi|jS/|rfw'|co2(K  -  Ko)G(0)(r',  r)G(0)(r',  r0) 

-  0~£)  [vlG(0,<r'’  r)'viG<oi(r'  r»)+s?G,o,(r'-  r)s?G,o)<r''  ■■»>]}• 


(A4) 


where  the  S'  integral  is  a  surface  integral  over  the  unper¬ 
turbed  interface  and  the  w’  integral  is  over  the  coordinate 
normal  to  the  unperturbed  interface.  The  range  of  the  w'  inte¬ 
gral  is  set  by  the  unperturbed  and  perturbed  interfaces  and  is 
[0,  £]  for  regions  where  the  perturbed  interface  lies  above 
the  unperturbed  interface  (£  >  0)  and  [£,  0]  for  the  other 
regions  (£  <  0).  In  the  integral  over  w',  the  zeroth-order 
Green’s  function  and  its  derivatives  are  approximated  as 
constant.  As  a  result,  the  w'  integral  introduces  an  overall 
factor  £,  from  which  it  follows  that  (A4)  gives  a  first-order 
correction  to  the  zeroth-order  field. 

In  Eq.  (A4),  the  transverse  and  normal  parts  of  the  gradi¬ 
ent  are  treated  separately  with  the  transverse  part  denoted 
V'j_ .  This  separation  is  employed  because  the  normal  deriva¬ 
tive  of  the  Green’s  function  is  discontinuous  at  the  unper¬ 
turbed  boundary.  In  contrast,  in  the  terms  not  involving 
normal  derivatives,  the  Green’s  function  and  its  transverse 
gradient  are  continuous  across  the  unperturbed  boundary  and 


may  be  taken  outside  the  W  integral  without  regard  to 
whether  they  are  evaluated  above  or  below  the  smooth 
interface. 

In  discussing  the  integral  over  w',  it  is  necessary  to 
break  the  integration  region  into  two  different  parts.  The 
“above”  regions  are  those  lens-shaped  regions  for  which 
£  >  0,  and  the  “below”  regions  are  those  for  which  £  <  0. 
Consistent  with  the  first-order  assumption,  all  acoustic 
parameters  are  taken  to  be  independent  of  w',  although  they 
may  be  functions  of  the  transverse  coordinates.  In  the 
“above”  regions,  k  —  k0  =  k2  —  Kq,  while  in  the  “below” 
regions,  k  —  k0  =  Ki  —  k2 .  Similar  expressions  hold  for  the 
density  factors.  As  noted  in  the  preceding  text,  the  integra¬ 
tion  limits  change  between  the  two  types  of  regions  with  the 
result  that  sign  differences  are  cancelled,  giving  factors  £(k2 
-  Ki )  and  £/(!/ p2  -  1/ pi)  for  both  regions. 

For  the  terms  in  Eq.  (A4)  that  do  not  involve  normal 
derivatives,  the  relevant  integrals  are,  therefore, 


FIG.  7.  Illustration  of  perturbation  of  smooth  bathymetry  by  a  small-scale 
roughness  contribution  measured  by  the  normal  displacement,  f(r'). 


(k  -  Ko)dw' =  {k2  -  Ki)C(r'),  (A5) 


and 


(A6) 


where  r'  is  evaluated  on  the  unperturbed  boundary,  S'. 
Because  the  integrals  are  proportional  to  £,  the  correspond¬ 
ing  terms  in  Eq.  (A4)  are  of  first  order,  as  anticipated. 

Because  of  the  discontinuity  at  w'  =  0  in  the  normal  de¬ 
rivative  of  the  zeroth-order  Green’s  function,  care  is  required 
in  treating  the  last  term  in  Eq.  (A4).  This  term  is  written 
in  an  ambiguous  fashion,  delaying  specification  of  whether 
the  normal  derivatives  are  to  be  evaluated  above  or  below 
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the  smooth  boundary.  The  first  normal-derivative  factor, 
r),  is  straightforward  and  is  evaluated  at  W 
=  0+  in  the  “above”  regions  and  at  w'  =  0“  in  the  “below” 
regions.  This  follows  simply  because,  in  the  starting 
Eq.  (A  1),  it  is  assumed  that  this  factor  is  to  be  evaluated  at 
the  integration  coordinate,  which  lies  above  the  smooth 
interface  in  the  “above”  regions  and  vice  versa.  In  the  first- 
order  approximation,  this  factor  is  assumed  constant  with 
respect  to  w\  but  it  takes  on  different  values  in  the  “above” 
and  “below”  regions. 

The  second  factor,  d/dw'G^°\ r',  r0)  is  more  subtle.  Two 
choices  are  possible,  each  leading  to  a  different  approxima¬ 
tion.  One  must  decide  whether  to  evaluate  this  derivative 
above  or  below  the  smooth  boundary.  It  should  be  remem¬ 
bered  that  one  seeks  an  approximation  to  the  derivative  of  the 
exact  Green’s  function.  For  the  “above”  regions,  extrapola¬ 
tion  from  below  should  result  in  a  reasonable  approximation 
to  the  derivative  over  the  region  of  integration,  as  this  region 
has  the  physical  properties  of  the  lower  medium.  For  the 
“below”  regions,  the  opposite  should  be  true.  Two  attractive 
properties  result  from  this  choice.  First,  reciprocity  is 
respected,  ps( r,  r0)/p(r)  =ps(r0,  r)/p(r0).  Second,  when  the 
zeroth-order  interface  is  flat,  the  result  is  found  to  be  identical 
to  that  given  by  conventional  small-roughness  perturbation 
theory. 

Defining 

'=J*'G“s)^G<0,(r'r)^G<0,(r'ro)' <A7> 


Eq.  (A  10).  The  source  position  r0  and  the  field  point  r  may  be 
separately  chosen  to  be  above  or  below  the  unperturbed 
boundary,  so  (4)  is  quite  general  and  can  be  applied,  e.g.,  to 
the  problem  of  scattering  from  the  water  into  the  water  as  well 
as  the  problem  of  scattering  from  the  water  into  the  sediment. 


APPENDIX  B:  PRE-SUMMED  ROUGH  SURFACE 
REALIZATIONS 

The  power  spectrum  for  the  2D  rough  surface  C(x,  y )  is 
assumed  to  be  of  the  form: 

IFC(K)  =  WC(KX1  Ky)  =  WC(K  cos  </>,  K sin</>),  (Bl) 

which  is  normalized  to  the  mean-square  height  of  the 
roughness  by  a2  =  J  d2KW^  (K) .  A  realization  of  the  two- 
dimensional  rough  surface  ((*,  y)  can  be  obtained  from  the 
two-dimensional  wave  vector  distribution  of  the  rough  sur¬ 
face  related  by  the  two-dimensional  Fourier  transform: 


C(R)  =  C(*,y) 

Z(K)  =  7T“T2 
{2n) 


=  d2kZ(  K)e,KR, 
J2RC(R)e“,K'R. 


(B2) 


The  values  of  the  distribution  at  two  different  wave  vector 
values  are  independent  Gaussian  random  variables,  and  the 
wave  vector  distributions  have  the  following  relation  to  the 
power  spectrum: 


one  finds  for  the  “above”  regions 


(Z*(K)Z(K'))  =  Wj(K)<5(K  —  K'), 


(B3) 


\p2  PiJ  dw 
X  ^7g(0)  (f/’  ro)|w'=°-C(r/), 


and  for  the  “below”  regions 


/=  f--1)  AG(o)(r>)|M/=0_ 

\Pi  PiJ  dW 
X  ^7^(0)(r,’  ro)U'=o+^(r,)• 


(A8) 


(A9) 


Expression  (A9)  can  be  brought  into  agreement  with 
Eq.  (A8)  through  use  of  the  displacement  matching  condition 


where  the  brackets  stand  for  ensemble  average.  Discretizing 
the  first  equation  of  (B2)  to  get  a  mean  value  of  the  distribu¬ 
tion  within  a  bin  of  wavenumbers: 


fR+()R 

z(K)  =  Z(K’)d2K', 

Jk 

and  using  Eq.  (B3) 


(B4) 


r 

(z*(K)z(K))  = 

Jk 

-f 


K+5K 


( Z*(K’)Z(K"))d2Kfd2K " 


K 
■K+<5K 


W{(K')8(K'  -  K  ")d2K'd2K" 


<Wc(K)AKxAKi 


■yi 


(B5) 


^Gl0|(r'r)l-'^  =^G<0)(r' r)'—-  <A10> 


a  realization  of  the  two-dimensional  surface  £(R)  can  be 
obtained  from  the  following  expression: 


The  transverse  and  normal  derivative  terms  can  be  combined 
into  a  single  three-dimensional  gradient  to  obtain  the  final 
result,  (4).  The  zeroth-order  Green’s  function  and  its  trans¬ 
verse  gradient  are  continuous  across  the  unperturbed  bound¬ 
ary.  Consequently,  they  can  be  evaluated  at  either  w  =  0+  or 
w  =  0~.  It  should  also  be  noted  that  the  choice  of  evaluating 
the  first  normal  derivative  factor  above  the  smooth  boundary 
and  the  second  below  is  arbitrary,  and  the  above  equation  is 
still  correct  if  the  opposite  choice  is  made.  This  follows  from 


(Nx-lNy-l  - 

C(mAx,  nAy)  =V2 Re  EEv  Wc(pAKx,qAKy)AKxAKy 

L  P=0  <7=0 

X  Rqp^PlN,+i2nnqlNy  \  ?  (B6 


where  m  =  0,1,2..., (Nx  -  1);  n  =  0,1,2..., (Ny  -  1)  and 
Rpq  =  [N(0, 1)  +  z7V(0, 1)]/V2,  with  N( 0,  1)  being  a  nor- 
mally  distributed  random  variable  with  zero  mean  and  unit 
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variance.  Multiplication  by  \fl  and  use  of  the  real  part  of  the  is  the  power  spectrum  averaged  over  the  azimuth.  Combin- 

complex  realization  yield  a  surface  having  the  desired  power  ing  Eqs.  (B9),  (BIO),  and  (B 11), 


spectrum  W^. 

Azimuthally  summed  one-dimensional  rough-surface 
realizations  can  be  generated  without  resorting  to  the  more 
numerically  demanding  generation  of  two-dimensional  real¬ 
izations.  As  noted  earlier,  the  azimuthally  summed  realiza¬ 
tions  are  defined  by 


n (R)  =  p£(fl,  <W, 

Jo 


(B7) 


and  it  is  required  that  the  equivalent  two-dimensional  real¬ 
izations  C(R,  (j))  have  the  specified  power  spectrum  W%( K). 
For  the  one-dimensional  surfaces,  there  exists  a  wavenumber 
distribution  through  Hankel  transforms 


i{R) 


p 

H(K)  = 

Jo 


KdKH(K)J0(KR), 


=  RdR  rj(R)J0(KR). 


(B8) 


This  distribution  is  also  statistically  independent  for  different 
wavenumbers: 


(H(K)H(K'))  =  Wn(K)S(K-K'). 


(B9) 


The  goal  now  is  find  a  relation  between  the  spectra  W^(K) 
and  W^K).  Using  Eq.  (B8),  the  following  relation  can  be 
obtained: 


RdR 


(H(K)H(K’))=  f 

J( 

x^RPiiR')). 

Using  Eqs.  (B7),  (B2),  and  (B3), 


R'dR'J0(KR)J0(K'R') 


(BIO) 


<; n{R)n{R ')>  = 


'2n  p2n 


d<t>  d(t>’ (£(x,y)^(x!  ,y')) 


0 

■2tt 


d(j)  d(j)r 

o  Jo  J 


dzK 


d2K'e~i{K  R- K' R,)  (Z\K)Z(K')) 


d1K 


p27 l 

2K'W((K)  5(K  —  K')  < 

Jo 


d(j) 


>2n 


d(j)  e 


i—iKR  cos {(j)-e)+iK'R'  cos {$ -&) 


(2k)2 1  d2KW/;( K)  J0(KR)J0(KR') 
(2k)3  | KdKWc( K)  J0(KR)J0(KR') 


(Bll) 


where 


■2n 


WC(K,  <j))d(j) 


(B12) 


Wt](K)8(K-K')  =  (H(K)H(K')) 


(2ti)3  Ji 


K"dK"Wc(K") 
RdRJ0(KR)J0(K"R) 
R'dR'  Jq  (K'R')Jq  (K"R') 
(2t r)3  \K"dK"W:(K") 


=  (2^) 


S(K"  -  K )  S(K"  -  K') 
K' 

8(K-K'). 


K 

3  WC(K) 

K 


(B13) 


Comparing  the  two  sides  of  the  preceding  equation,  the 
desired  relation  between  W^{k)  and  Wv(k)  is  obtained: 


Wn(K)  =  (2n) 


,WC(K) 


K 


(B14) 


To  generate  numerical  realizations  of  the  azimuthally 
summed  random  surface,  the  second  equation  in  Eq.  (B8) 
must  be  converted  to  discrete  form.  To  that  end,  define  the 
following  discrete  distribution  within  a  bin  of  wavenumber: 


K+AK 


h(K)  = 
From  Eq.  (B9), 

(h2(K))  = 


dK'H(K'). 


(B15) 


?K+AK 

pK+AK 

dK' 

K 

K 

?K+AK 

pK+AK 

dK' 

K 

K 

dK"  (H(K')H(K")) 
\iK"W,i(K')5(K’-K") 


■  W,I(K)AK. 


Let 


h(K)  =  ^JwJk)AK  N(0, 1). 


(B16) 


(B17) 


Then  a  realization  can  be  generated  by  the  following  digi¬ 
tized  version  of  the  first  equation  of  (B8): 


»/C R )  =Y,h(Kn)KnMK„R), 


(B18) 


where  Kn  =  nAK,  n  =  0,  1,2,  ....  Let  the  maximum  range  of 
the  surface  be  L/2,  and  A KL  =  2n, 


N- 1 


r]  ( mAR)=J2h(Kn)KnJ0  ( KnmAR ) 

n= 0 

=  (2rc)2El 


n=0 


jWr(nAK)(nAK 

^ - An 

V  L  1 

iAKAR)Rt 
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where  Rin  =  A(0,  1),  and  Eq.  (B14)  has  been  used  to  express 
the  realization  using  To  estimate  the  power  spectrum 
from  realizations  of  rough  surfaces,  use: 


Wn(K)  = 


A  K 


{H2(K))AK 


=  2i(H\K)) 


2n 


POO 

RdR  rj(R)Jo(KR) 

Jo 


(B20) 


APPENDIX  C:  CALCULATING  BACKSCATTERING 
STRENGTH 

When  an  ensemble  of  backscattered  time  series  ps(t)  is 
obtained  either  from  measurement  or  from  simulation,  the  back- 
scattering  cross  section  is  calculated  according  to  the  definition: 


(Cl) 


where  I0(t)  is  the  incident  intensity  at  a  unit  distance  from 
the  source  and  dA  =  2nRdR  is  the  ensonified  area,  where  R 
is  the  horizontal  range.  The  same  expression  also  applies  for 
the  slant  range,  r,  dA  =  2nrdr.  It  is  assumed  that  the  source 
has  no  azimuthal  directivity.  The  average  denoted  by  ()  is 


over  a  finite  sample  of  realizations,  unlike  the  infinite  ensem¬ 
ble  employed  in  Appendix  B.  Because  r  =  c\t/ 2,  dr  =  c\dt/ 2 
and  J  I0(t)dt  =  E  is  the  source  energy, 


(Woi2)(ci03 

%tic\E 


(C  2) 


APPENDIX  D:  BACKSCATTERING  STRENGTH  IN 
PERTURBATION  THEORY 

It  will  be  shown  that  eq.  (10)  leads  to  a  backscattering 
cross  section  expression  that  is  the  same  as  given  by  first- 
order  perturbation  theory.  First,  rewrite  Eq.  (10)  in  its  two- 
dimensional  form 


k2  f°° 

p-{z°)=iL 


„2ik\r' 


dx' dy' q(x'  ,/)  — —  t  (Rf).  (Dl) 


Assume  that  the  rough  patch  of  area  A  is  centered  at  (xc,  yc). 
Then 

k\r  «  k\rc  +  SR  •  <5K,  (D2) 

where  rc  =  y/x2  +  y2c  +zg,  <5R  =  (x'  -  xc,yf  -  yc)  =  R'  -  Rc, 
<5K  =  (k\  cos0ccos</>c,  k\  cos0c  sin(/)c),  sin0c  =  zo/rc,  and 
cos  (j)c  =  (xf  —  xc)/\5R\.  Then  the  mean  scatterered  intensity 
is  approximately 


to)r 


~  (27T)2 

=  %V{Rc)\2dA 


POO  P( 

dx’dy’ 

— oo  J— oo  J 


2ik\{r'-r") 

dx"dy"(ax',y'W,  /)) - 2—  |T(i?c)  |2 


{2n) 


—  OO 
'OO  POO 


POO 

d{x'  -  x")d{y'  -  y")N{x'  -  x" ,  /  - 

—  OO  J— OO 


=  k2\x(Rc)\2dAW(2dK), 


(D3) 


where  N  is  the  two-dimensional  correlation  function  and  W  its 
power  spectmm.  According  to  the  definition  of  Appendix  C, 
the  backscattering  cross  section  is 

°(0)  =  (\Ps(zo)\2)r2c/dA 

=  ki\T(Rc)\2W(2kicosOccos(j)c:2ki  cos0csin(/>c),  (D4) 
in  agreement  with  the  usual  expression21  (Chap.  13). 
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Geoacoustic  inversion  work  has  typically  been  carried  out  at  frequencies  below  1  kHz,  assuming 
flat,  horizontally  stratified  bottom  models.  Despite  the  relevance  to  Navy  sonar  systems  many  of 
which  operate  at  mid-frequencies  (1-10  kHz),  limited  inversion  work  has  been  carried  out  in  this 
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bottom  loss  data  between  2  and  5  kHz.  The  acoustic  measurements  were  taken  during  the  Shallow 
Water  2006  Experiment  off  the  coast  of  New  Jersey.  A  half-space  bottom  model,  with  three 
parameters  density,  compressional  wave  speed,  and  attenuation,  was  used  for  inversion  by  fitting 
the  model  to  data  in  the  least- square  sense.  Inverted  sediment  sound  speed  and  attenuation  were 
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ing  the  known  spatial  variability  in  this  area.  The  observations  and  modeling  results  demonstrate 
that  forward  scattering  from  topographical  changes  is  important  at  mid-frequencies  and  should  be 
taken  into  account  in  sound  propagation  predictions  and  geoacoustic  inversion.  To  cope  with 
fine-scale  topographic  variability,  measurement  technique  such  as  averaging  over  tracks  may  be 
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I.  INTRODUCTION 

Sediment  geoacoustic  parameters,  in  many  cases,  are 
the  most  important  parameters  for  predicting  shallow  water 
sound  propagation  and  reverberation.  Due  to  the  difficulty  of 
in  situ  measurements,  indirect  methods  have  been  commonly 
used  to  invert  for  sediment  geoacoustic  parameters.  Among 
various  geoacoustic  inversion  techniques,  bottom  reflection/ 
loss  measurements  have  been  used  to  extract  sediment  geoa¬ 
coustic  properties. 

Since  the  1960  s,  studies  on  bottom  reflection  have  been 
carried  out1-3  to  understand  the  interaction  between  sound 
and  the  ocean  bottom.  Comparisons  have  been  made 
between  measured  values  of  reflection  coefficients  (or  bot¬ 
tom  loss)  and  plane  wave  reflection  coefficients  with  the  bot¬ 
tom  modeled  as  a  layered-fluid  medium.1-6  A  greater  part  of 
these  studies  were  at  low  frequencies,  i.e.,  below  1  kHz,1,4-6 
with  considerable  acoustic  energy  penetrating  the  bottom. 
This  results  in  additional  complexity  in  interpreting  bottom 
reflection  measurements  if  bottom  properties  vary  signifi¬ 
cantly  with  depth.  Such  is  the  case  with  a  so-called  transition 
layer  where  bottom  properties  vary  linearly  with  depth.  Due 
to  the  sound  speed  gradient,  energy  refracted  within  the  bot¬ 
tom  can  be  significant  in  modeling  bottom  loss  (BL)  data  at 
low  frequencies.4-6  Analytic  solutions  have  been  derived  for 
the  reflection  of  plane  acoustic  waves  from  horizontally 
stratified  fluid  layer  whose  density  and  sound  speed  vary 
linearly  with  depth.7-9 

As  to  geoacoustic  inversion  using  bottom  reflection, 
some  of  the  published  work  focused  on  estimating  sediment 
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properties  within  a  transition  layer,  i.e.,  to  invert  for  gra¬ 
dients  of  sound  speed  or  attenuation  within  the  layer.  Using 
a  fluid  bottom,  Spofford10  proposed  a  technique  to  estimate 
the  sound  speed  gradient  using  acoustic  data  between 
50-1600  Hz.  Holland  et  al.n  demonstrated  a  method  to 
obtain  density  and  sound  speed  gradients  in  the  transition 
layer  (top  1.5  m)  by  a  Bayesian  inversion  approach  in  the 
frequency  band  of  300-1600  Hz.  Many  inversion  models 
assume  discrete  layering.  A  unique  inversion  scheme  was 
presented  by  Holland  et  al.12  based  on  a  joint  time-  and 
frequency-domain  data  analysis  using  wideband  signals 
between  600  and  6000  Hz.  The  time-domain  BL  data  were 
used  to  obtain  layer  thicknesses  and  interval  velocities, 
whereas  the  frequency-domain  data  were  used  to  determine 
depth-  and  frequency-dependent  attenuation  and  depth- 
dependent  density.  In  addition,  geoacoustic  inversion  using 
BL  data  was  studied  where  the  sediment  was  modeled  as  an 
elastic13  or  poro-elastic14,15  medium.  Since  a  fluid  bottom  is 
adopted  in  the  present  work,  the  details  of  the  elastic  and 
poro-elastic  models  will  not  be  addressed  here. 

As  noted  above,  geoacoustic  inversion  work  has  typi¬ 
cally  been  carried  out  at  low  frequencies.  Despite  the  rele¬ 
vance  to  Navy  sonar  systems,  many  of  which  operate  at  mid¬ 
frequencies  (1-10  kHz),  limited  inversion  work  has  been  car¬ 
ried  out  in  this  frequency  band  where  sound  propagation  and 
reverberation  can  be  strongly  influenced  by  variability  in  the 
sea  bottom,  surface,  and  water  column.  This  paper  is  an 
effort  to  demonstrate  the  viability  of  geoacoustic  inversion 
using  bottom  loss  data  at  mid-frequencies  (2-5  kHz).  The 
acoustic  data  presented  in  this  work  were  taken  during  the 
Shallow  Water  2006  Experiment  (SW06)  off  the  coast  of 
New  Jersey.16  As  one  of  the  mid- frequency  efforts  in  SW06, 
the  acoustic  measurement  was  conducted  in  the  central  area 
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A.  Tow  data 


FIG.  1.  (Color  online)  GPS  locations  of  acoustic  measurements.  Diamond: 
moored  receiving  array;  black  line:  acoustic  track  with  a  towed  source. 


of  SW06  (location  shown  in  Fig.  1).  During  SW06,  many 
acoustic  measurements  were  carried  out  in  this  area,  and  one 
of  the  main  goals  is  to  improve  inversion  technique  by  com¬ 
paring  geoacoustic  inversion  results  using  different  methods 
with  supporting  ground  truth  measurements.  The  ground 
truth  measurements  include  several  earlier  experiments  and 
a  geological  survey  conducted  on  mid  and  outer  New  Jersey 
shelf  to  characterize  the  seabed17-21  and  in  situ  measure¬ 
ments  in  SW06. 22-24  We  demonstrate  that  it  is  possible  to 
use  mid-frequency  data  to  successfully  estimate  sediment 
sound  speed;  however,  it  is  also  found  that  care  needs  to  be 
taken  to  minimize  the  impact  of  fine- scale  topographical 
variation. 

This  paper  is  organized  as  follows.  In  Sec.  II,  the  acous¬ 
tic  data  and  results  are  presented,  and  the  bottom  modeling 
and  geoacoustic  inversions  results  are  shown  in  Sec.  III. 
Section  IV  summarizes  the  main  results  and  discusses  the 
implications  of  this  work. 

II.  ACOUSTIC  DATA  AND  RESULTS 

The  acoustic  data  involved  in  this  analysis  were 
obtained  along  a  track  with  a  slowly  towed  source.  Data 
were  taken  on  August  19,  2006  during  the  SW06  experiment 
and  recorded  on  a  receiving  array  moored  at  39°  01.47'  N, 
73°  02.26'  W.  On  the  receiving  array,  there  are  two  sub¬ 
arrays,  clustered  at  depths  25  and  50  m  below  the  surface, 
and  each  has  four  elements.  For  the  sub-array  at  25  m,  the 
four  elements  are  at  25.0,  25.2,  25.5,  and  26.4  m;  for  the  one 
at  50  m,  the  four  elements  are  at  50.0,  50.2,  50.5,  and  51.4  m. 
Linear- frequency-modulated  (LFM)  chirp  signals  in  the  fre¬ 
quency  band  of  1.5-10.5  kHz  were  transmitted  using  an 
ITC-2015  acoustic  source  (International  Transducer  Com¬ 
pany)  with  10%  cosine  taper  at  the  beginning  and  end  of  the 
chirp.  The  focus  of  this  work  is  on  data  received  at  25  m. 


The  acoustic  data  were  taken  on  August  19, 
00:29:00-00:49:00  UTC.  The  source  was  towed  by  the  R /V 
Knorr  at  a  speed  of  0.5  knot  along  an  80  m  isobath,  and  LFM 
signals  were  transmitted  every  20  s  with  a  source  depth  of 
40  m.  The  tow  started  at  a  range  of  105  m  and  ended  at 
approximately  8  km.  The  data  shown  in  this  work  correspond 
to  the  beginning  part  of  the  tow,  i.e.,  from  105  m  to  400  m 
(Fig.  1).  There  are  a  total  of  60  pings  during  this  tow 
interval. 

The  signals  received  at  25  m  were  first  pulse- 
compressed  in  the  frequency  band  of  1.5-6  kHz.  In  Fig.  2, 
the  compressed  signals  for  the  25  m  receiver  are  plotted 
against  the  ping  number  on  the  vertical  axis  and  the  reduced 
travel  time  on  the  horizontal  axis  with  a  time  equal  to  propa¬ 
gation  time  for  1495  m/s  at  the  range  between  the  source  and 
receiving  array  for  each  ping  removed.  The  first  three  arriv¬ 
als  were  identified  as  direct  arrival,  surface  and  bottom 
reflections.  Since  our  goal  is  to  study  bottom  loss,  the  focus 
will  be  on  the  bottom  reflection.  As  shown  in  Fig.  2(b),  all 
bottom  bounces  were  aligned  at  time  zero  and  cut  with  a 
2.5  ms  time  window.  The  time  series  were  carefully  exam¬ 
ined  and  there  was  no  evidence  of  secondary  reflectors  in  the 
bottom  reflection  signal.  The  60  pings  shown  here  corre¬ 
spond  to  grazing  angles  from  10°  to  approximately  43°. 

To  obtain  BL,  a  sensitivity  calibration  is  needed,  and  the 
source  beam  pattern  has  to  be  taken  into  account.  Here,  the 
direct  path  of  the  tow  data  will  be  used  for  calibration.  The 
advantage  of  using  the  direct  path  can  be  seen  from  Fig.  3 
which  shows  the  sound  speed  profile  on  the  left  and  geome¬ 
try  of  acoustic  measurements  on  the  right.  When  towed  at  a 
constant  speed,  the  source  was  tilted  and  this  tilt  angle,  0,  is 
the  same  for  both  the  direct  path  and  the  bottom  bounce 
path.  This  will  simplify  beam  pattern  correction  due  to  the 
tilt.  In  Fig.  3,  r0  is  the  distance  between  the  ship  stern  and 
the  receiving  array  obtained  from  GPS  and  ship  gyro  record¬ 
ings;  v i  is  the  propagation  distance  for  the  bottom  bounce, 
obtained  from  arrival  time.  Using  the  geometry  shown  in 
Fig.  3,  the  tilt  angle  can  be  computed  via 


10  20  30  40 

Reduced  travel  time  (ms) 


FIG.  2.  (Color  online)  (a)  Pulse  compressed  results  for  60  pings  with 
reduced  travel  time  (see  definition  in  text)  in  the  frequency  band  of 
1.5-6  kHz;  (b)  same  data  with  the  bottom  bounce  lined  up  at  time  zero. 
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f  2H  —Z,=ZS  cos  9  +  \fr\~-X 
\  y  +  Zs  sin  9  =  r0 


(1) 


where  Zr  is  the  receiver  depth,  Zs  is  the  distance  from  the 
ship  stern  to  the  source,  and  H  is  the  water  depth.  It  is 
assumed  that  the  tow  line  is  straight,  that  is,  that  the  effect  of 
drag  on  the  tow  line  is  negligible  compared  to  the  effect  of 
the  tension  force  exerted  by  the  source.  The  tilt  angle 
obtained  from  Eq.  (1)  for  each  ping  falls  in  the  angle  range 
of  1.5°  ±  0.5°  throughout  the  acoustic  measurements.  Conse¬ 
quently,  a  constant  value  of  1.5°  was  applied  throughout  the 
tow. 

As  mentioned,  direct  path  data  will  be  used  to  calibrate 
the  bottom  bounce.  Similarly  to  the  bottom  bounces,  the  first 
19  direct-path  signals  were  lined  up  and  cut  with  a  2.5  ms  time 
window.  Figure  4  shows  their  receive  levels  which  were  cor¬ 
rected  for  spreading  and  the  source  beam  pattern  (see  Sec.  II  B) 
at  2,  3,  4,  and  5  kHz.  Transmissions  corresponding  to  the  short¬ 
est  ranges  (at  highest  grazing  angles)  appear  on  the  right-hand 
side  of  the  plot.  The  receive  levels  of  the  first  7  pings,  corre¬ 
sponding  to  a  grazing  angle  range  of  10.7°-9.6°,  are  relatively 
stable  with  variation  ±0.5  dB.  For  pings  further  down  range, 
i.e.,  at  smaller  grazing  angles,  the  intensity  variation  can  be  up 
to  5  dB.  This  is  due  to  the  fact  that  the  direct  path  straddles  the 
thermocline  and  slight  changes  in  the  water  column  can  have 
significant  effects  on  this  path  even  at  a  propagation  distance  of 
approximately  150  m  (Ref.  25).  The  first  7  pings  were  used  to 
obtain  an  average  to  calibrate  all  bottom  bounce  signals.  The 
average  signal  was  obtained  by  first  lining  up  the  7  pings  and 
taking  a  coherent  average. 

B.  Source  beam  pattern 

To  obtain  BL,  it  is  critical  to  correct  both  the  direct  path 
and  the  bottom  bounce  path  for  the  source  beam  pattern.  Fac¬ 
tors  that  affect  the  beam  pattern  correction,  such  as  the  source 
tilt  and  azimuthal  rotation,  have  been  carefully  examined. 
Three-dimensional  beam  pattern  measurements  at  2,  3,  4,  and 
5  kHz  were  carried  out  at  the  APF  dockside  calibration  facil¬ 
ity.  Results  shown  in  Fig.  5  are  vertical  beam  patterns  at  one 


FIG.  4.  (Color  online)  Direct  path  intensity  corrected  for  spreading  and  the 
source  beam  pattern  at  2,  3,  4,  and  5  kHz  (pings  1-19  have  decreasing  graz¬ 
ing  angles). 

of  the  several  azimuth  angles  in  the  measurement  set  with  the 
source  sitting  upright.  Specifically,  the  figure  shows  the  beam 
pattern  at  the  azimuth  angle  which  faces  the  receiver  under 
hydrodynamic  considerations  of  the  experimental  setup.  The 
direct  path  from  the  source  to  the  receiver  at  25  m  depth  cor¬ 
responds  to  the  angular  range  of  272°-281°  in  the  beam  pat¬ 
terns.  Figure  5  shows  considerable  anti-symmetry  in  the 
vertical,  i.e.,  top  versus  bottom,  due  to  a  user-installed  suspen¬ 
sion  system  at  the  bottom  of  the  source.  The  3-D  beam  pattern 
data  will  be  used  in  identifying  the  beam  pattern  uncertainty 
due  to  azimuthal  rotation  of  the  source. 


C.  Bottom  loss 

Bottom  loss  can  be  obtained  via 


BL  =  10  x  log10 


'WW/flfo  -  e)\ 

„  tiixdSiixd/Dtii  +  0)  /  ’ 


(2) 


where  I0(Xo)  and  h( 7j)  are  intensities  of  the  calibration  and 
the  ith  ping;  S0(Xo)  and  S[( xd  are  corresponding  spreading 
losses  obtained  through  ray  tracing;  and  D(xo~S)  and 


Ship 


FIG.  3.  (Color  online)  (left)  Sound  speed  profile  for  the  tow;  (right)  geometry  for  the  towed  source  data. 

J.  Acoust.  Soc.  Am.,  Vol.  131 ,  No.  2,  Pt.  2,  February  2012  Yang  etal Mid-frequency  geoacoustic  inversion  1713 


Downloaded  27  Feb  2012  to  128.95.76.55.  Redistribution  subject  to  ASA  license  or  copyright;  see  http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp 


2  kHz 
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0 


FIG.  5.  (Color  online)  Measured  source 
beam  pattern  at  2,  3,4,  and  5  kHz  at  one  azi¬ 
muth  angle  with  source  sitting  upright.  The 
direct  path  from  the  source  to  the  receiver  at 
25  m  depth  corresponds  to  the  angular  range 
of  272°-281°. 


D(/j  +  8)  are  source  beam  pattern  corrections  with  tilt  angle 
taken  into  account.  Both  y0  and  yr  are  launch  angles  for  cali¬ 
bration  and  BL  measurement,  respectively. 

BL  results  at  four  frequencies,  2,  3,  4,  and  5  kHz,  are 
shown  in  Fig.  6  for  a  10%  bandwidth.  Results  at  different 
frequencies  are  quite  close  to  each  other,  and  the  critical 
angle  is  well  defined  (around  25°  for  all  frequencies).  The 
result  at  2  kHz  for  grazing  angles  above  35°  is  about  2dB 
lower  than  the  other  three,  which  is  due  to  insufficient  water 
depth  in  the  beam  pattern  measurements  leading  to  a  prob- 
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lem  with  multipathing.  In  addition,  all  four  curves  show  a 
feature  around  20°.  Modeling  work  was  carried  out  to  under¬ 
stand  the  feature  and  possible  explanations  involve  local  var¬ 
iations  in  (1)  sediment  properties;  (2)  bottom  layering;  and  (3) 
topography.  The  modeling  results  will  be  shown  in  Sec.  Ill 
together  with  the  inversion  results. 


D.  Uncertainties  in  BL  data 

Consider  a  general  case  with  measurement  variable  Y 
where  Y  is  a  function  of  j  variables: 


Y  =  Y(XuX2,...,Xj). 

Then  the  uncertainty  of  measurement  variable  Y  can  be  writ¬ 
ten  as26 


UY= 


2 


M~2UX2)  +"'  + 


1/2 


(3) 


where  Uxt  is  the  uncertainty  for  the  zth  variable,  and  where 
the  errors  in  each  variable  are  assumed  to  be  independent. 
Following  (3),  we  can  write  the  uncertainty  for  BL: 


£4  =  K  +  ujo  +  u2m)  +  uj  +  u2m .  (4) 


FIG.  6.  (Color  online)  Bottom  loss  results  at  2,  3,  4,  and  5  kHz  with  a  10% 
bandwidth. 


UfQ  is  the  uncertainty  of  the  direct  path  intensity.  Recall  that 
the  direct  paths  of  the  first  seven  pings  have  been  used  to 


1714  J.  Acoust.  Soc.  Am.,  Vol.  131 ,  No.  2,  Pt.  2,  February  2012 


Yang  etal.\  Mid-frequency  geoacoustic  inversion 


Downloaded  27  Feb  2012  to  128.95.76.55.  Redistribution  subject  to  ASA  license  or  copyright;  see  http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp 


2  kHz 


3  kHz 


4  kHz 


5  kHz 


Grazing  angle  (°) 


FIG.  7.  (Color  online)  Bottom  loss  results  with 
uncertainty  bounds  at  2,  3,  4,  and  5  kHz. 


find  an  average  signal  for  calibration,  therefore,  UjQ  is  taken 
as  the  variance  of  the  intensity  of  the  first  seven  pings. 

UjQ  and  U\  are  uncertainties  of  the  spreading  loss  for 
the  direct  and  bottom-bounce  paths  modeled  by  ray  tracing 
due  to  limited  knowledge  of  the  environment.  Since  the 
source  was  towed  and  the  receiving  array  did  not  have 
oceanographic  sensors,  an  oceanographic  mooring,  approxi¬ 
mately  1.5  km  away  from  the  receiving  array,  was  used  to 
provide  sound  speed  profiles  for  modeling.  There  are  10  sen¬ 
sors  on  the  mooring  and  the  sampling  interval  is  30  s  (sensor 
depths:  13.3,  15.1,  18.8,  22.6,  26.3,  33.8,  41.3,  56.3,  71.3, 

78.5  m).  The  60  pings  were  recorded  between  00:29:00  and 
00:49:00  UTC  on  August  19.  Within  the  20-min  recording 
time,  there  were  a  total  of  39  profiles  recorded  by  the  ocean¬ 
ographic  mooring.  Since  the  environmental  measurement 
was  not  done  along  the  acoustic  track,  the  question  now  is, 
how  much  error  is  incurred  in  modeling  spreading  loss  if  dif¬ 
ferent  sound  speed  profiles  are  used?  Individual  profiles  are 
used  to  model  the  intensity  of  both  the  direct  paths  and  bot¬ 
tom  bounces.  At  each  frequency,  six  cases  have  been  tried 
using  individual  profiles  No.  5,  10,  15,  20,  30,  and  a  20-min 
average  profile.  UjQ  and  Uj  are  defined  as  the  standard  devi¬ 
ations  of  the  results  (not  shown  here)  from  these  six  cases. 

For  beam  pattern  uncertainties,  U and  U^ey  each 
have  three  components.  The  three  components  are  source  tilt 
angle  uncertainty  due  to  variations  in  acoustically  deter¬ 
mined  values,  tilt  angle  uncertainty  due  to  GPS  resolution, 
and  source  azimuthal-rotation-induced  uncertainty.  The  tilt 
angle  computed  using  propagation  distances  varies  between 

1.5  ±  0.5°  and  the  0.5°  uncertainty  in  tilt  angle  was  trans¬ 
lated  into  uncertainty  using  beam  pattern  data.  As  shown 
earlier,  GPS  recordings  were  used  to  find  distances  between 
the  source  and  receiver  to  compute  the  tilt  angle.  In  this 
experiment,  GPS  has  a  ±  0.5  m  resolution  which  adds  extra 

J.  Acoust.  Soc.  Am.,  Vol.  131 ,  No.  2,  Pt.  2,  February  2012 


tilt  angle  uncertainty  which  can  be  similarly  converted  into 
uncertainty  using  beam  pattern  data.  For  source  rotation 
induced  uncertainty,  it  was  assumed  that  the  source  can 
rotate  ±90°.  Three-dimensional  beam  pattern  data  were 
used  to  provide  uncertainty  due  to  the  rotation. 

With  all  the  uncertainties  identified,  the  BL  results  are 
re-plotted  in  Fig.  7  at  2,  3,  4,  and  5  kHz  with  uncertainty 
bars.  The  uncertainties  at  3  and  4  kHz  are  smaller  than  those 
at  2  and  5  kHz,  which  is  due  to  the  fact  that  uncertainty  in 
beam  pattern  correction  dominates.  In  addition,  the  uncer¬ 
tainty  bound  itself  is  smaller  than  the  BL  fluctuations.  We 
consider  these  fluctuations  may  be  real  and  possibly  due  to 
forward  scattering  from  topographical  variation.27  To  take 
out  the  fluctuations,  one  would  need  to  repeat  acoustic  meas¬ 
urements  over  multiple  tracks  and  do  ensemble  averaging. 

III.  BOTTOM  MODELING  AND  GEOACOUSTIC 
INVERSION 

A  simple  half  space  bottom  model,  with  three  parame¬ 
ters,  density,  compressional  wave  speed,  and  attenuation, 
was  used  for  geoacoustic  inversion  by  best  matching  model 
with  data  in  the  least  square  sense.  The  choice  of  a  half¬ 
space  over  a  layered  bottom  model,  which  was  also  investi¬ 
gated  but  not  shown  here,  is  based  on  the  fact  that  little 
frequency  dependence  was  observed  in  BL  measurements 
(Fig.  6).  One  expects  frequency-dependent  interference  with 
a  layered  bottom,  therefore,  a  half-space  model  is  appropri¬ 
ate  for  the  present  work.  In  addition,  BL  obtained  with  a  flat, 
plane  wave  model  was  used  to  compare  with  data  for  inver¬ 
sion.  Considering  the  receiver  depth  above  the  water/sedi¬ 
ment  interface  (55  m)  relative  to  the  wavelength  (0.75  m), 
the  approximation  of  geometrical  acoustics  is  satisfied  and 
the  plane-wave  reflection  coefficient  can  be  used.28  Inverted 

Yang  etal .:  Mid-frequency  geoacoustic  inversion  1715 


Downloaded  27  Feb  2012  to  128.95.76.55.  Redistribution  subject  to  ASA  license  or  copyright;  see  http://asadl.org/journals/doc/ASALIB-home/info/terms.jsp 


E 

o 

O) 

CL 


1500  1600  1700  1  2  3  4  1  2  3  4 

cb  (m/s)  ab  (dB/m)  ab  (dB/m) 


1 

0.9 

0.8 

0.7 

0.6 

0.5 

0.4 

0.3 

0.2 

0.1 


FIG.  8.  Cost  function  surfaces  in  terms  of  p-c,  c-a,  and  p-ot  through  an  exhaustive  search  at  3  kHz. 


sediment  geoacoustic  parameters  and  their  uncertainties  are 
summarized  in  Sec.  Ill  A.  Additional  modeling  work  has 
been  carried  out  in  an  attempt  to  understand  the  interesting 
feature  around  20°  in  BL.  Possible  explanations  for  the  fea¬ 
ture  and  modeling  results  are  presented  in  Sec.  Ill  B.  Direct 
measurement  should  be  used  as  ground  truth  to  validate 
inversion  results.  Therefore,  in  Sec.  Ill  C,  inverted  sediment 
sound  speed  results  of  this  work  are  compared  with  direct 
measurements  and  other  inversion  results  which  were  carried 
out  within  a  2  km  x  2  km  boxed  area. 

A.  Half-space  fluid  model 

A  simple  half  space  bottom  model  is  used  with  three  pa¬ 
rameters,  density  p ,  compressional  wave  speed  c,  and  attenu¬ 
ation  a,  and  they  are  allowed  to  be  frequency  dependent.  An 
exhaustive  search  was  carried  out  to  estimate  the  three  pa¬ 
rameters  by  best  matching  model  with  data  in  the  least 
square  sense.  One  advantage  of  doing  an  exhaustive  search 
is  that  we  can  look  at  correlations  between  different  parame¬ 
ters.  The  cost  function  is  defined  as  the  squared  error 
between  model  and  data  summed  over  all  angles.  Figure  8 
shows  the  cost  function  surfaces  in  terms  of  p-c,  c-a,  and  p- oc 
at  3  kHz,  and  the  red  area  indicates  the  most  probable  region 
for  the  inversion  results.  The  correlation  between  parameters 
can  be  clearly  seen  from  Fig.  8:  p-c  has  an  inverse  relation 
since  the  product  of  the  two  determines  acoustic  impedance; 
c  and  oc  are  slightly  correlated;29  while  p  and  oc  are  quite 


independent  of  each  other,  i.e.,  there  is  little  correlation 
between  the  two. 

The  best  fits  between  data  and  a  half-space  model  at  2, 
3,4,  and  5  kHz  are  shown  in  Fig.  9,  and  the  inversion  results 
for  p,  c,  and  oc  and  their  uncertainties  are  summarized  in  Fig.  10. 
The  uncertainties  were  computed  as  follows.30  For  a  linear  sys¬ 
tem  of  equations  d  =  Gm,  the  least  squares  solution  is 

mL2  =  (CfGf'c/d. 


For  a  weakly  nonlinear  problem  like  the  one  in  this  paper, 
we  can  expand  the  problem  around  the  solution  point  to  the 
first  order,  and  therefore,  matrix  G  can  be  defined  as 
Gp  =  Specifically,  for  BL  with  three  parameters  p,  c, 


and  oc,  i.e.,  m  = 


P 

c 

oc 


,  and  N  data  points,  G  can  be  written  as 


dBL\  (m) 

0BLt(m) 

dBL\  (m) 

dp 

dc 

da 

dBL2(m ) 

dBL2(m ) 

dBL2(m) 

dp 

dc 

doc 

dBLN(m ) 

dBLN(m ) 

dBLN(m) 

dp 

dc 

doc 

where  ^  =  201og10(e)Re{/T  ||}  ^  and  R  is  the  complex 
reflection  coefficient. 
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FIG.  9.  (Color  online)  Best  match 
between  model  and  data  at  (a)  2  kHz, 
(b)  3  kHz,  (c)  4  kHz,  (d)  5  kHz  includ¬ 
ing  the  feature  around  20°. 


The  system  of  equations  is  scaled  by  W  =  diag(l/cri, 
1/(72,  •••  1/ Gj)  and  the  least  squares  solution  is  then 
mLi  =  (GTwGw)~lGTwdw,  where  Gw  =  WG  and  dw  =  Wd.  The 
covariance  of  a  data  vector  d  of  normally  distributed,  independ¬ 
ent  random  variables  with  linear  transformation  matrix  A  is 

Co  \(Ad)  =  ACov(d)AT. 

With  A  =  (GtwGw)  '  Gtw,  the  covariance  matrix  for  mLl  can 
be  written  as 

Cov(mL2)  =  (GtwGw)  1GtwIGw  (GtwGw)~X 

=  {GtwGw)~\  (6) 

The  uncertainty  is  defined  as  the  95%  confidence  interval  for 
each  model  parameter  and  can  be  obtained  via  the  covari¬ 
ance  matrix: 

mLi±\.96  ■  diag(Cov(mL2))1/2.  (7) 

Using  (7),  the  uncertainties  for  p,  c,  and  a  can  be  obtained. 
Figure  10(a)  summarizes  the  inversion  results  using  the 
entire  BL  data  and  Fig.  10(b)  using  the  BL  data  excluding 
the  feature  around  20°.  Figure  10  shows  that  in  general, 
uncertainty  increases  with  frequency  for  all  frequencies.  For 
density  and  sound  speed,  the  uncertainties  are  quite  small 
while  for  attenuation,  the  uncertainty  can  sometimes  be  up 
to  50%.  By  comparing  Figs.  10(a)  and  10(b),  it  is  clear  that 
inverted  density  and  sound  speed  change  little  with  or  with¬ 
out  the  feature  around  20°  while  attenuation  has  been  greatly 
reduced.  Based  on  the  work  of  Williams  et  al .,31  1%  disper¬ 
sion  in  sound  speed  is  expected  in  the  frequency  band  of 
2-5  kHz.  Dispersion  at  this  level  would  be  masked  by  uncer¬ 
tainty  in  sediment  sound  speed  for  this  work. 

J.  Acoust.  Soc.  Am.,  Vol.  131 ,  No.  2,  Pt.  2,  February  2012 


As  known,  BL  has  different  sensitivity  to  p,  c,  and  a  at 
different  grazing  angles.  For  grazing  angles  above  the  criti¬ 
cal,  BL  is  sensitive  to  density  while  for  those  below  the  criti¬ 
cal,  BL  is  sensitive  to  attenuation.  The  critical  angle  depends 
on  sound  speed;  therefore  BL  is  sensitive  to  sound  speed  for 
angles  near  critical.  As  a  result,  one  expects  little  change  in 
inversion  results  for  p  and  c  with  or  without  the  feature 
around  20°  while  the  contrary  is  expected  for  attenuation. 

B.  Possible  explanation  for  feature  around  20° 

As  mentioned  in  Sec.  II,  possible  explanations  for  the 
feature  around  20°  are  local  variations  in  (1)  sediment  prop¬ 
erty;  (2)  interference  due  to  a  layered  bottom;  and  (3)  topog¬ 
raphy.  For  possibility  (1),  a  local  soft  spot  is  needed  to 
reproduce  the  feature.  Based  on  a  geological  survey  on  the 
New  Jersey  middle  and  outer  shelf,18  this  explanation  is 
unlikely.  The  98  grab  samples  and  collocated  acoustic  meas¬ 
urements  using  the  In  situ  Sound  Speed  and  Attenuation 
Probe  (ISSAP)  show  that  the  sediment  at  the  SW06  experi¬ 
mental  site  is  comprised  mostly  of  clayey  sands  to  coarse 
sands.  Specifically  for  the  location  of  the  measurements  in 
this  work,  the  sediment  is  classified  as  “ribbons”18  with  a 
mean  grain  size  between  1  (/>  and  2.5  0  with  no  examples  of 
the  very  fine-grained  sediment,  i.e.,  grain  size  around  4  </>, 
required  to  explain  the  feature.32  For  possibility  (2),  a  model 
with  a  layer  of  finite  thickness  over  a  semi-infinite  layer  was 
used  to  invert  for  local  parameters.  The  model  has  seven  pa¬ 
rameters:  2  sets  of  p,  c,  and  a  for  the  finite  and  semi-infinite 
layers  respectively  and  a  layer  thickness.  Inversion  results 
(not  shown  here)  yield  layer  thickness  that  is  inversely  pro¬ 
portional  to  frequency,  ruling  out  the  hypothesis  of  local  var¬ 
iation  in  bottom  layering. 

To  investigate  possibility  (3),  a  broadband  simulation 
was  carried  out  to  model  the  experiment  using  the  parabolic 
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FIG.  10.  (Color  online)  Summary  of  inversion  results  for  p,  c,  and  a  (a)  using  entire  BL  data  and  (b)  using  BL  data  excluding  feature  around  20°. 


equation  (PE)  with  bathymetry  variation.  To  reproduce  the 
feature,  this  bathymetry  variation  was  found  to  span  about 
50  m  resembling  a  convex  surface,  and  a  Gaussian  shape  is 
used  to  represent  this  bathymetry  change.  This  feature  was 
positioned  to  cover  the  range  interval  of  90-140  m  from  the 
source  corresponding  to  grazing  angles  between  16°-24°. 
Figure  11(a)  shows  the  simulation  results  for  BL  with  a  flat 
bottom  for  reference  and  Fig.  11(b)  with  an  8-cm-high  Gaus¬ 
sian  bottom  feature.  The  simulation  results  show  a  BL  fea¬ 
ture  similar  to  that  observed  in  data,  and  the  width  and 
magnitude  match  the  experiment  quite  well.  The  simulation 
results  are  surprising  in  amplitude  considering  such  a  gentle 
topographic  change.  A  simulation  using  the  Kirchhoff 


approximation  was  used  to  check  the  results  obtained  from 
PE.  The  results  from  PE  and  Kirchhoff  are  in  good  agree¬ 
ment,  demonstrating  a  maximum  2  dB  BL  change  due  to  the 
bottom  feature  with  reference  to  a  flat  bottom  case.  Both  PE 
and  Kirchhoff  are  2D  simulations  and  no  out-of-plane  scat¬ 
tering  is  allowed.  In  addition,  simulation  results  show  that 
this  bottom  feature  affects  acoustics  mainly  at  low  grazing 
angles. 

Based  on  the  modeling  results,  possibility  (3)  may  be 
the  cause  for  the  feature  around  20°.  The  simulation  results, 
however,  are  unable  to  reproduce  the  slight  frequency  de¬ 
pendence  of  the  feature  observed  in  data  (Fig.  6)  due  to  the 
fact  that  it  was  hypothesized  as  a  geometrical  lensing  effect. 


FIG.  11.  (Color  online)  PE  modeling  results  at  2,  3,  4,  and  5  kHz  with  (a)  flat  bottom  and  (b)  an  8  cm  high  Gaussian  bottom  feature. 
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FIG.  12.  Comparison  of  sediment  sound  speed  between  direct  measurements  and  geoacoustic  inversion  results  in  SW06.  (a)  GPS  locations  for  each  measure¬ 
ment  and  (b)  sediment  sound  speed  results.  Water  sound  speed  near  water/sediment  interface  is  approximately  1498  m/s  for  SW06. 


This  frequency  dependence  may  require  further  modeling 
effort  in  the  future  and  can  be  due  to  finer- scale  topographi¬ 
cal  changes  than  the  Gaussian  feature.  Unfortunately,  there 
is  no  ship  echo  sounder  data  or  bathymetry  data  at  the  mea¬ 
surement  site  with  enough  resolution  to  confirm  the  topo¬ 
graphical  change.  However,  a  bottom  feature  such  as  that 
shown  here  could  occur  even  for  a  relatively  flat  site  such  as 
this.  The  importance  of  the  modeling  work  is  that  it  indicates 
how  detailed  the  environmental  sampling  should  be  for  mid¬ 
frequency  propagation  and  inversion  studies.  In  cases  where 
it  is  not  practical  to  use  detailed  bathymetry  in  inversion, 
averaging  over  several  tracks  may  compensate  the  error  seen 
in  the  present  work.27 

C.  Comparison  of  sediment  sound  speed  and 
attenuation  between  direct  measurement  and 
inversion 

Direct  measurement  and  geoacoustic  inversion  results  of 
sediment  sound  speed  and  attenuation  in  SW06  are  com¬ 
pared  within  a  2  km  boxed  area.  The  GPS  plot  for  all  meas¬ 
urements  is  shown  in  Fig.  12(a)  and  the  sediment  sound 
speed  results  are  shown  in  Fig.  12(b).  The  sediment  attenua¬ 
tion  results  are  summarized  in  Table  I. 

In  the  2  km  boxed  area,  three  direct  measurements  and 
three  geoacoustic  inversions  are  included.  The  three  direct 
measurements  use  different  instruments  and  they  are  Sedi¬ 
ment  Acoustic- speed  Measurement  System  (SAMS,  Ref.  24), 
ISSAP22,23  and  Geo  Probe.33  The  penetration  depths  of 
SAMS,  ISSAP,  and  Geo  Probe  are  1.6,  0.2,  and  0.3  m, 
and  measurements  were  carried  out  at  2-20,  65,  and  20  kHz. 
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The  three  inversion  results  included  are  the  one  from  this 
work,  work  from  Jiang  et  al?4  and  Choi  et  al.?5  respec¬ 
tively.  The  three  inversion  results  are  from  short  range  prop¬ 
agation  data  in  the  frequency  bands  of  2-5,  1.5-4. 5,  and 
1-20  kHz,  respectively. 

For  sediment  sound  speed  comparison  (Fig.  12)  of  the 
three  direct  measurement  results,  the  ones  using  SAMS  and 
Geo  Probe  are  within  each  other’s  uncertainty  bounds.  Two 
of  the  inversion  results,  Jiang  et  al 34  and  Choi  et  al.35  are 
consistent  with  these  two  direct  measurements.  In  addition, 
all  four  measurements  as  shown  in  Fig.  12(a)  are  in  close  vi¬ 
cinity.  The  other  two  results,  i.e.,  the  direct  measurement 
result  using  ISSAP  and  the  inversion  result  of  this  work, 
give  higher  sound  speed  and  yet  do  not  seem  mutually  con¬ 
sistent  given  the  uncertainty  in  each.  Note  that  the  ISSAP 
operates  at  65  kHz  and,  accounting  for  the  3%  dispersion31 


TABLE  I.  Comparison  of  sediment  attenuation  between  direct  measure¬ 
ment  and  geoacoustic  inversion  within  a  2  km  boxed  area  in  SW06. 


Reference  name 

Operating  Frequency 
(kHz) 

Attenuation 

(dB/m/kHz) 

ISSAP 

65 

0.60 

Geo  probe 

20 

0.22 

Yang,  Jackson, 
and  Tang 

2,  3,  4,  and  5 

[0.05  0.17  0.45  0.48] 
(BL  exclude  feature) 
[0.20  0.40  0.63  0.60] 
(entire  BL) 

Jiang,  Chapman, 
and  Gerstoft 

1.75-3.15 

0.1 

Choi,  Dahl,  and  Goff 

1-20 

0.05 
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expected  between  65  and  3  kHz,  the  IS  SAP  result  would  be 
consistent  with  the  present  work.  The  discrepancy  between 
four  measurements  with  lower  values  and  the  other  two,  con¬ 
sidering  their  geographical  locations,  is  likely  due  to  spatial 
variation  of  sediment  properties.  Substantial  spatial  variation 
in  sediment  sound  speed  over  small  spatial  scales  has  been 

reported  in  this  area  through  the  ONR-sponsored  Geoclutter 
1 8 

program. 

Sediment  attenuation  results  from  the  same  measure¬ 
ments  in  Fig.  12  except  SAMS  are  summarized  in  Table  I.  In 
Table  I,  attenuation  results  of  this  work  at  individual  fre¬ 
quencies,  i.e.,  2,  3,  4,  and  5  kHz  are  listed  for  two  cases:  esti¬ 
mations  using  BL  data  with  and  without  the  20°  feature. 
Table  I  shows  that  attenuation  of  five  measurements  ranges 
from  0.05-0.6  dB/m/kHz  within  the  2  km  boxed  area.  It  is  of 
interest  to  compare  results  from  this  work  with  those  from 
Jiang  et  al.  and  Choi  et  al.  which  are  in  the  same  frequency 
range.  The  attenuation  results  of  this  work  are  consistent 
with  Jiang  et  al.  and  Choi  et  al.  if  BL  data  excluding  the  20° 
feature  are  used  for  inversion.  As  mentioned  earlier,  BL  is 
sensitive  to  attenuation  at  low  grazing  angles,  where  BL  data 
are  more  susceptible  to  topographical  changes.  Therefore, 
the  estimation  of  attenuation  is  more  affected  by  forward 
scattering  from  topographical  changes  than  sediment  density 
and  sound  speed  (Fig.  10).  For  fine-scale  topographic  vari¬ 
ability,  the  estimation  of  attenuation  can  be  improved  by 
repeating  the  experiment  along  different  tracks  and  carrying 
out  ensemble  averaging. 

IV.  SUMMARY  AND  DISCUSSION 

In  this  paper,  geoacoustic  inversion  results  using  bottom 
loss  data  collected  from  SW06  in  the  frequency  range  of 
2-5  kHz  have  been  presented.  Measured  BL  data  exhibited  a 
well-defined  critical  angle  at  all  frequencies  and  little  fre¬ 
quency  dependence  was  observed.  Uncertainties  in  BL  were 
carefully  identified,  and  it  was  found  that  uncertainties  due 
to  beam  patterns,  excluding  fine-scale  topography,  are  the 
dominant  component. 

A  simple  half  space  fluid  model,  with  three  parameters 
p,  c,  and  a,  was  used  for  geoacoustic  inversion.  Inversion 
results  were  obtained  by  best  matching  between  model  and 
data  in  the  least  square  sense.  Inverted  sediment  sound 
speeds  were  used  to  compare  with  direct  measurements  and 
results  using  other  inversion  techniques  in  a  2  km  x  2  km 
area  in  SW06.  The  comparison  yields  consistency  and  dis¬ 
crepancy  which  are  both  correlated  with  the  geographical 
locations  for  the  measurements.  This  leads  to  the  conclusion 
that  the  observed  discrepancy  may  be  due  to  spatial  varia¬ 
tions  of  sediment  properties.  The  comparison  of  attenuation 
obtained  in  the  same  frequency  range  show  that  they  are  con¬ 
sistent  and  within  the  uncertainty  bounds  of  this  work  if  the 
feature  at  20°  is  excluded  for  inversion. 

The  present  work  demonstrates  the  practicality  of 
inverting  geoacoustic  parameters  using  BL  data  at  mid¬ 
frequencies.  In  addition,  the  observations  and  modeling 
results  from  this  work  show  that  forward  scattering  from  to¬ 
pographical  changes  is  important  at  mid-frequencies  and 
should  be  taken  into  account  in  sound  propagation  predic- 
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tions  and  geoacoustic  inversion.  Most  inversion  work 
assumes  flat,  horizontally  stratified  bottom  models  and  scat¬ 
tering  from  the  bottom  is  usually  ignored.  In  this  work,  the 
fluctuations  observed  in  BL  exceed  the  uncertainty  bounds 
and  may  be  due  to  forward  scattering  from  fine-scale  topo¬ 
graphical  changes.  The  feature  around  20°  may  also  be  due 
to  topographical  change.  Modeling  work  carried  out  to 
explain  the  feature  requires  a  gentle  but  consistent  topo¬ 
graphical  change  with  elevation  on  the  order  of  10  cm  and 
spatial  span  of  50  m.  Though  this  is  a  small  elevation  change, 
the  resultant  BL  variation  is  about  2  dB  with  reference  to  a 
flat  surface  case  between  2  and  5  kHz.  In  the  present  work, 
fluctuations  in  BL  due  to  forward  scattering  from  topograph¬ 
ical  changes  do  not  alter  the  inversion  results  for  sediment 
sound  speed  and  density  significantly  but  do  affect  the  deter¬ 
mination  of  attenuation  and  uncertainty  bounds  for  all 
inverted  parameters.  In  addition  to  uncertainty  analysis  of 
geoacoustic  inversion,  the  results  of  this  work  also  have 
impact  on  the  development  of  measurement  methods  that 
can  cope  with  fine-scale  topographic  variability,  e.g.,  averag¬ 
ing  over  tracks. 
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Abstract:  Knowledge  of  sediment  sound  speed  is  crucial  for  predicting 
sound  propagation.  During  the  Shallow  Water  ’06  experiment,  in  situ  sedi¬ 
ment  sound  speed  was  measured  using  the  Sediment  Acoustic-speed  Mea¬ 
surement  System  (SAMS).  SAMS  consists  of  ten  fixed  sources  and  one  re¬ 
ceiver  that  can  reach  a  maximal  sediment  depth  of  3  m.  Measurements  were 
made  in  the  frequency  range  2-35  kHz.  Signal  arrival  times  and  propagation 
distances  were  recorded,  from  which  sediment  sound  speed  was  determined. 
Preliminary  results  from  three  deployments  show  that  SAMS  was  capable  of 
determining  sediment  sound  speed  with  uncertainties  less  than  1.6%.  Little 
dispersion  in  sediment  sound  speed  was  observed. 
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1.  Introduction 

Acoustic  interaction  with  the  sea  bottom  is,  in  many  cases,  an  essential  component  of  sound 
propagation  in  a  shallow  water  waveguide.  In  situ  geoacoustic  properties  of  the  seabed,  how¬ 
ever,  are  difficult  to  obtain.  Since  the  1950s,  in  situ  direct  measurements  have  been  carried  out 
to  study  the  sediment  geotechnical  properties.1  More  recently,  sediment  geoacoustic  properties 
have  been  measured  within  the  surfacial  layer  in  the  frequency  range  1-100  kHz  using  either 
manually  buried  acoustic  systems  (Refs.  2-4)  or  specially  designed  underwater  mechanical 
systems.  Refs.  5-8,  from  the  latter  category,  present  systems  that  can  penetrate  into  the  sedi¬ 
ment  through  their  own  gravitational  forces.  The  In  Situ  Sediment  geoacoustic  Measurement 
System  (Refs.  5-7)  is  designed  for  in  situ  measurements  within  the  topmost  30  cm  while  the 
acoustic  lance  (Ref.  8)  has  a  maximum  penetration  depth  of  5  m. 

As  part  of  the  experimental  effort  in  Shallow  Water  ’06  (SW06),  the  Sediment 
Acoustic-speed  Measurement  System  (SAMS)  was  used  to  directly  measure  the  sediment 
sound  speed.  SAMS  is  driven  into  the  seabed  by  a  powerful  vibrocore,  which  allows  precise 
penetration  depth  up  to  3  m  with  arbitrary  step  size.  The  ground  truth  measurements  are  valu¬ 
able  not  only  in  studying  in  situ  sediment  properties  but  also  in  providing  sediment  geoacoustic 
data  to  which  inversion  results  can  be  compared. 

This  paper  is  organized  as  follows.  In  Sec.  2,  the  analyses  of  calibration  and  sediment 
data  are  presented.  Section  3  summarizes  and  discusses  future  directions.  System  uncertainty 
analysis  is  given  in  the  Appendix. 

2.  Data  analysis 

SAMS  was  deployed  during  the  SW06  field  experiment.  Four  data  sets,  one  in  the  water  column 
as  calibration  and  three  in  the  sediment,  were  recorded.  For  all  data  sets,  three  linear- frequency- 
modulated  (LFM)  “chirps”  were  used  in  the  frequency  bands  2-11,  10-21,  and  20-35  kHz, 
which  are  referred  to  as  low-,  mid-,  and  high-frequency  (LF,  MF,  and  HF)  in  later  analysis. 

2.1  Analysis  of  calibration  data 

Calibration  data  were  acquired  in  a  bay  environment  and  relatively  far  from  the  positions  of  the 
sediment  measurements.  For  calibration,  SAMS  was  deployed  such  that  both  sources  and  re- 
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Time  (ms)  Residual  (|as) 

Fig.  1.  (Color  online)  Pulse  compression  results  at  MF.  (a)  Waveforms  received  from  source  10;  (b)  linear  regression 
results;  (c)  linear  regression  residual  scatter  plot;  (d)  residual  histogram. 


ceiver  were  always  within  the  water  column.  The  analysis  of  the  calibration  data  is  used  to 
augment  the  system  geometry  measurements  and  establish  system  uncertainty. 

SAMS  consists  of  a  3-m-tall  triangular  frame  and  a  2-m-long  extension  beam.  Three 
sources  are  fixed  on  the  triangular  frame  with  the  other  seven  on  the  extension  beam.  The  re¬ 
ceiver  is  located  in  the  center  of  the  frame  and  can  be  driven  vertically  into  the  sediment  using 
a  vibrocore.  Horizontal  distances  between  sources  and  receiver  are  between  1.18  and  2.97  m. 

Data  were  taken  at  3 1  receiver  depths  with  a  stepsize  of  about  0.1  m.  At  each  depth,  the 
ten  sources  sequentially  transmitted  three  chirps  and  each  was  repeated  five  times.  Pulse  com¬ 
pression  was  carried  out  as  an  initial  step  in  processing  the  data.  In  Fig.  1(a),  pulse  compressed 
waveforms  received  from  source  ten  are  plotted  at  their  corresponding  receiver  depths.  Signal 
arrival  time  is  defined  as  the  time  at  the  peak  of  the  envelope  and  highlighted  with  a  dot.  With 
signal  arrival  times  defined,  the  speed  of  sound  in  water  can  be  determined  from  the  linear 
regression  of  the  arrival  times  and  distances  between  the  sources  and  receiver.  Five  repeated 
pings  are  used  to  find  the  averaged  arrival  times  yielding  a  total  of  310  data  points.  The  linear 
regression  result,  Fig.  1(b),  is  1502.2  m/s  with  2.7  m/s  uncertainty  at  the  95%  confidence 
level.  The  calculation  of  uncertainty  assumes  that  the  residuals  (difference  between  data  and  fit 
function)  are  random  and  follow  a  normal  distribution  of  zero  mean  and  constant  variance.  A 
scatter  plot  of  the  residuals  is  shown  in  Fig.  1(c).  Both  the  scatter  plot  and  its  histogram  [Fig. 
1(d)]  indicate  that  the  distribution  of  the  residual  is  close  to  normal.  Similar  procedures  are 
repeated  for  the  LF  and  HF  calibration  data.  The  curve  fitting  results  for  the  speed  of  sound  in 
water  are  1503. 4 ±6.7  and  1503.1  ±3.4  m / s,  respectively. 

The  sound  speeds  determined  are  close  to  each  other  with  the  confidence  interval  at 
MF  and  HF  almost  completely  enclosed  by  that  of  the  LF.  The  higher  uncertainty  at  LF  is  due  to 
a  roughly  58%  decrease  in  the  total  number  of  data  points  included  in  the  curve  fitting  process. 
Beyond  2  m,  the  LF  calibration  data  showed  signals,  possibly  due  to  the  tube  waves,  arriving 
prior  to  the  direct  arrivals  and,  therefore,  were  excluded  from  the  curve  fitting.  Unfortunately, 
there  was  no  conductivity-temperature-depth  (CTD)  record  when  the  calibration  data  were 
taken.  Concurrent  ship  data  only  provide  temperature  and  salinity  at  the  sea  surface.  Therefore, 
a  total  of  120  historical  summer  CTD  data  around  the  area  were  sought  for  reference.  Individual 
distance  to  the  SAMS  calibration  position  varied  from  16  to  52  km  and  water  depth  ranged 
13-40  m.  A  strong  thermocline  was  observed  starting  at  around  10  m  and  water  sound  speed 
varied  from  1520  m/s  at  the  sea  surface  to  1495  m/s  at  35  m.  For  calibration,  bathymetry 
showed  a  23  m  water  depth  with  SAMS  suspended  5-6  m  below  the  surface,  i.e.,  measure¬ 
ments  were  taken  at  6-9  m  in  depth.  Compared  with  historical  CTD  data,  the  calibration  re¬ 
sults,  1503.4±6.7,  1502.2±2.7,  and  1503.1  ±3.4  m/s,  respectively,  are  reasonably  within  the 
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Fig.  2.  (Color  online)  MF  Pulse-compressed  waveforms  from:  (a)  source  1;  (b)  source  10  at  different  receiver  depths, 
(c)  and  (d)  are  half-space  Green’s  function  simulations  of  (a)  and  (b),  respectively. 


variation  range.  More  importantly,  the  acoustic  data  are  of  good  quality  and  repeatable  at  dif¬ 
ferent  frequency  bands,  supporting  the  conclusion  that  the  water  sound  speeds  measured  by 
SAMS  are  close  to  the  true  values  at  the  time. 

2.2  Analysis  of  sediment  data 

To  start,  two  cases  in  the  MF  band  from  position  2  are  chosen  to  show  the  characteristics  of 
sediment  data.  The  first  case,  Fig.  2(a),  draws  the  signals  from  source  1  at  19  receiver  depths 
with  the  maximum  penetration  depth  around  1.6  m.  Data  are  processed  in  the  same  way  as  the 
calibration.  The  receiver  starts  at  about  1 0  cm  above  the  sediment  surface,  which  makes  the  first 
depth  sample  less  than  zero  in  the  figure.  In  determining  the  signal  arrival  time,  a  time  window 
is  specified  (bounded  by  the  black  dots)  within  which  the  peak  of  the  signal  envelope  is  recog¬ 
nized  and  highlighted  with  a  green  dot.  It  is  quite  obvious  that  there  are  two  erroneous  readings 
of  the  peak  time  in  Fig.  2(a)  as  the  receiver  first  enters  the  sediment.  For  sources  that  are  further 
away  from  the  receiver,  there  are  more  such  occurrences  as  shown  in  Fig.  2(b).  Data  of  this  kind 
are  carefully  excluded  from  analysis.  Signal-to-noise  ratio  drops  considerably  for  the  geom¬ 
etries  realized  in  Fig.  2(b)  due  to  the  combined  effects  of  longer  acoustic  path  through  the 
sediment  and  ray  bending  by  the  critical  angle.  A  two-half  space  Green’s  function  is  used  to 
simulate  the  scenario  of  Figs.  2(a)  and  2(b)  and  the  results  are  shown  in  Figs.  2(c)  and  2(d).  The 
critical  angle  effect  is  quite  apparent  in  Fig.  2(d). 


EL118  J.  Acoust.  Soc.  Am.,  Vol.  124,  No.  3,  Pt.  2,  September  2008  Yang  et  al.:  Direct  measurement  of  sediment  sound  speed 


Downloaded  07  Oct  2013  to  128.95.76.55.  Redistribution  subject  to  ASA  license  or  copyright;  see  http://asadl.org/terms 


Yang  et  al:  JASA  Express  Letters 


[DOI:  10.1121/1.2963038] 


Published  Online  28  August  2008 


Table  1.  Summary  of  sediment  sound  speed  results. 


Cb  (m/s) 

Position  1 

LF 

1614.8  ±8.7 

MF 

1622.1  ±12.5 

Position  2 

LF 

1597.7  ±11.0 

MF 

1598.6  ±9.8 

Position  3 

LF 

1588.2  ±  15.8 

MF 

1611. 6±  24.8 

To  find  sediment  sound  speed,  the  water  sound  speed  close  to  the  water-sediment  in¬ 
terface  is  required.  Throughout  the  experiment,  CTD  records  show  a  very  stable  water  sound 
speed  at  depths  beyond  70  m.  Based  on  CTD  records,  the  water  sound  speed  is  1496  m/s. 
Sediment  sound  speed  is  assumed  homogeneous  within  the  penetrated  depth,  1.6  m.  Ray  trac¬ 
ing  is  carried  out  for  each  data  point  by  varying  sediment  sound  speed  in  the  model  and  the 
closest  match  of  arrival  times  between  measurement  and  ray  tracing  determines  the  in  situ  sedi¬ 
ment  sound  speed.  The  uncertainties  are  determined  in  the  same  manner  as  for  the  calibration 
(Sec.  2.1).  Assuming  uncertainty  comes  entirely  from  propagation  in  the  sediment,  the  time  and 
distance  that  are  spent  in  water  are  removed  from  the  total  time  and  distance.  The  uncertainty  in 
sediment  sound  speed  is  then  calculated  using  the  residuals  at  95%  confidence  level. 

Results  are  summarized  in  Table  1  and  Fig.  3.  Results  from  positions  1  and  2  show 
similar  uncertainty  bounds  around  10  m/s  and  little  dispersion  between  the  two  frequency 
bands.  There  is  about  a  20-m/ s  sound  speed  difference  between  the  two  positions,  which  is 
believed  to  be  the  true  spatial  variation.  Results  at  position  3  have  much  larger  uncertainties 
than  the  other  two.  Signs  of  signal  degradation  were  observed  (not  shown  here),  which  results  in 
a  50%  higher  mean  residual  than  the  other  two  positions.  Geological  studies  around  the  SW06 
region  (Ref.  9),  using  interpreted  chirp  seismic  reflection  data,  indicate  a  substantial  difference 
in  sediment  properties  at  position  3  from  positions  1  and  2,  which  is  believed  to  be  the  cause  for 
higher  uncertainty  at  position  3.  In  addition,  acoustic  measurements  of  bottom  reflection  (Ref. 
1 0)  were  made  in  the  vicinity  of  SAMS  positions  1  and  2  in  the  frequency  range  1-20  kHz.  The 
geoacoustic  inversion  results  of  sediment  sound  speed,  with  co-located  coring  and  stratigraphic 
studies,  are  consistent  with  the  direct  measurement  results  using  SAMS. 

3.  Summary  and  future  directions 

In  this  paper,  direct  measurements  of  sediment  sound  speed  using  the  Sediment  Acoustic-speed 
Measurement  System  (SAMS)  have  been  presented.  The  calibration  data  were  first  analyzed  to 
establish  the  system  uncertainty,  which  is  approximately  3  m/s.  Sediment  data  were  taken  at 
three  positions.  Sediment  sound  speeds  and  uncertainties  are  summarized  in  Table  1  and  Fig.  3. 
Results  indicate  a  20-m/s  sound  speed  variation  between  positions  1  and  2.  At  position  3,  the 


Fig.  3.  (Color  online)  Sediment  sound  speeds  with  uncertainties  at  three  positions. 
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increase  in  uncertainty  may  be  attributed  to  the  sediment  properties  based  on  geological  studies 
around  the  central  experimental  area  in  SW06.  The  sediment  sound  speeds  found  at  positions  1, 
2,  and  3arel618±ll,1598±10,  and  1600  ±20  m/s,  respectively  Little  dispersion  in  sediment 
sound  speed  was  observed.  Direct  measurement  of  sediment  sound  speed  dispersion  has  been 
found  in  sandy  sediments  (Refs.  2  and  11).  The  dispersion  was  observed  to  be  at  its  greatest  in 
the  frequency  range  800-2000  Hz.  In  this  work,  the  frequency  coverage  is  higher  than  the 
pronounced  transition  region  of  dispersion,  which  may  explain  the  observed  lack  of  dispersion. 
Future  directions  include  improvement  of  system  uncertainty  in  sediment  sound  speed,  deter¬ 
mination  of  sediment  attenuation  and  its  dispersion  relation,  and  depth  dependence  of  sediment 
geoacoustic  properties. 
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Appendix:  System  uncertainty  analysis 

For  a  system  like  the  SAMS,  the  dimension  and  propagation  time  uncertainties  limit  the  reso¬ 
lution  of  sound  speed  measurement.  Uncertainty  can  come  from  both  measurements  and  meth¬ 
odology  utilized  to  analyze  data.  In  this  work,  part  of  the  uncertainty  comes  from  measurement 
of  distances,  i.e.,  horizontal  distance  and  initial  depth  offset  between  sources  and  receiver;  ini¬ 
tial  depth  offset  between  receiver  and  the  sediment  surface;  and  receiver  depth  reading.  The 
reading  of  the  arrival  times  falls  in  the  latter  category. 

Following  Ref.  12,  a  general  function  q  with  multiple  variables  (x,y,z,...)  has  uncer¬ 
tainty: 


In  Eq.  (Al),  variables  x,y,z,...  are  independent  measurements  with  uncertainties 
Sx,  Sy,  Sz, . . ..  For  this  work,  the  general  function  is  the  speed  of  sound  in  water,  cw,  which  is  the 
ratio  of  distance  and  time: 


r(x,y,z,  •  •  • )  Vx2  +  (h  +  d )2 


(A2) 


where  r  is  the  slant  distance,  t  is  travel  time;  x,  h,  and  d  are  horizontal  distance,  receiver  depth, 
and  initial  vertical  distance  between  source  and  receiver.  Following  Eq.  (Al),  the  uncertainty  in 
calibration  can  be  written  as 

■  (a3) 

The  individual  uncertainties  Sx,Sh ,  Sd  are  defined  as  0.9,  0. 1 ,  and  0.5  cm.  Specifically,  Sx  is  set 
to  \  of  the  source  dimension;  Sd  is  set  to  \  of  the  receiver  dimension;  Sh  accounts  for  depth 
reading  uncertainty.  In  Eq.  (A3),  the  most  difficult  part  is  to  determine  St.  The  receiver  is  em¬ 
bedded  inside  a  stainless  steel  tube  with  two  rectangular  windows  open  on  the  side.  The  com¬ 
bination  of  direct  arrival  and  reflections  off  the  window  may  slightly  change  the  location  of  the 
signal  peak.  Assuming  this  window  effect  is  random,  St  can  be  determined  using  the  mean 
residual  obtained  from  Fig.  1  as  7.5  jus.  Figure  4  shows  the  uncertainty  corresponding  to  each 
source.  For  each  of  them,  as  depth  increases,  the  uncertainty  decreases  from  top  to  bottom.  The 
black  dashed  line  is  the  mean  value  for  each  source.  It  is  obvious  that  sources  that  are  closer  to 
the  receiver  have  higher  uncertainty  and  spreading.  The  least  uncertainty  for  an  individual  mea¬ 
surement  is  around  5  m/s. 

In  calibration  analysis,  data  recorded  from  all  ten  sources  were  used  in  the  linear  re¬ 
gression  to  find  the  speed  of  sound  in  water,  i.e.,  a  relationship  between  individual  (Fig.  4)  and 
overall  uncertainty  (Fig.  1)  needs  to  be  clarified.  The  overall  system  uncertainty  can  be  deter- 
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Fig.  4.  (Color  online)  System  uncertainty  in  the  speed  of  sound  in  water. 


mined  from  the  individual  uncertainties  in  Fig.  4  as  a  two-step  process.  First,  assume  an  ideal 
system  with  no  uncertainty,  i.e.,  find  propagation  times  by  dividing  measured  propagation  dis¬ 
tances  by  a  fixed  sound  speed.  The  value  cw=  1502.2  m/s,  determined  from  calibration  at  MF,  is 
used.  The  linear  regression  shows  a  perfect  fit  between  propagation  distance  and  time  with  zero 
uncertainty.  Second,  convert  the  maximum  uncertainty  of  each  source  Scw ,  as  in  Fig.  4,  to  its 
equivalent  uncertainty  in  distance,  i.e.,  by  multiplying  Scw  with  corresponding  propagation 
time  t.  Then,  the  measured  propagation  distances  were  added  a  random  quantity  in  the  range  of 
±tX  Scw.  The  linear  regression  is  carried  out  again  by  forcing  cw=  1502.2  m/s.  The  system 
uncertainty  at  95%  confidence  interval  is  calculated  in  a  similar  fashion  as  in  Sec.  2.  This  un¬ 
certainty,  determined  to  be  ±1.6  m/s,  is  the  overall  system  uncertainty. 
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Effect  of  the  Internal  Tide  on  Acoustic  Transmission 

Loss  at  Midfrequencies 
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Abstract — Nonlinear  internal  waves  are  a  common  event  on  the 
continental  shelf.  The  waves  depress  the  high-gradient  region  of  the 
thermocline  and  thicken  the  surface  mixed  layer  with  consequent 
effect  on  acoustic  propagation.  After  the  waves  have  passed,  it  may 
take  several  hours  for  the  thermocline  to  rise  to  its  prewave  level. 
To  examine  the  effect  of  the  rising  thermocline,  oceanographic  and 
acoustic  data  collected  during  the  2006  Shallow  Water  Experiment 
(SW06)  are  analyzed.  Midfrequency  acoustic  data  (1.5-10.5  kHz) 
taken  for  several  hours  at  both  fixed  range  (550  m)  and  along  a  tow 
track  (0.1-8.1  km)  are  studied.  At  the  fixed  range,  the  rising  ther¬ 
mocline  is  shown  to  increase  acoustic  intensity  by  approximately 
5  dB.  Along  the  tow  track,  the  transmission  loss  changes  2  dB  for  a 
source-receiver  pair  that  straddles  the  thermocline.  Using  oceano¬ 
graphic  moorings  up  to  2.2  km  away  from  the  acoustic  receiver,  a 
model  for  the  rising  thermocline  is  developed.  This  ocean  model 
is  used  as  input  to  a  broadband  acoustic  model.  Results  from  the 
combined  model  are  shown  to  be  in  good  agreement  with  experi¬ 
mental  observation.  The  effects  on  acoustic  signals  are  shown  to  be 
observable,  significant,  and  predictable. 

Index  Terms — Acoustic  signal  processing,  nonlinear  ocean  in¬ 
ternal  waves,  underwater  acoustic  telemetry. 


I.  Introduction 

DENSITY  stratification  in  the  ocean  leads  to  the  propaga¬ 
tion  of  internal  waves.  In  deep  water,  these  waves  are  es¬ 
sentially  linear  and  are  commonly  modeled  as  a  random  process 
[1].  In  shallow  water,  more  event-like,  nonlinear  waves  [2]  can 
be  generated  that  display  bore-like  features  [3] .  Driven  by  the  in¬ 
teraction  of  the  tide  with  bottom  topography,  nonlinear  internal 
waves  depress  the  high-gradient  region  of  the  thermocline  on 
the  order  of  10  m  and  propagate  at  wave  speeds  on  the  order 
of  1  m/s.  The  effects  these  propagating  ocean  waves  have  on 
acoustic  propagation  can  be  dramatic.  At  low  frequencies,  they 
produce  an  observed  rapid  decorrelation  of  acoustic  modes  [4] 
and  a  subsequent  recorrelation  that  can  be  explained  using  a 
coupled  mode  model  [5],  [6].  Also  at  low  frequencies,  they  can 
cause  horizontal  refraction  leading  to  focusing  and  defocusing 
[7],  [8].  The  presence  of  nonlinear  internal  waves  affects  the 
horizontal  coherence  length  at  both  low  [9]  and  high  [10],  [11] 
frequencies. 
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In  the  cited  papers,  the  primary  emphasis  was  on  how  acoustic 
signals  are  affected  by  nearby  nonlinear  internal  waves.  The  in¬ 
ternal  waves,  however,  are  but  one  manifestation  of  the  internal 
tide  that  is  generated  by  the  tide  passing  over  the  continental 
shelf  break.  Due  to  nonlinearities,  the  deeper  thermocline  part  of 
the  internal  tide  travels  faster  than  the  higher  thermocline  part. 
In  a  rather  short  distance,  the  deepest  part  of  the  wave  catches  up 
to  the  shallowest  part  and  high-frequency  nonlinear  waves  form 
as  part  of  the  resulting  bore.  A  second  manifestation  of  the  in¬ 
ternal  tide  is  the  slow  rising  of  the  thermocline  after  the  wave 
has  passed.  In  this  paper,  the  emphasis  is  on  how  acoustic  sig¬ 
nals  are  affected  by  these  gradual  changes  in  the  water  column. 

The  topic  is  relevant  to  acoustics  because  the  rising  thermo¬ 
cline  caused  by  the  internal  tide  may  last  for  several  hours  while 
nonlinear  internal  waves  may  transect  a  given  acoustic  track  for 
only  a  small  fraction  of  the  tidal  period.  At  short  ranges,  where 
the  acoustic  arrivals  can  be  separated  in  time,  the  gradual  rising 
of  the  thermocline  is  shown  to  change  acoustic  intensity  by  ap¬ 
proximately  5  dB.  At  longer  ranges,  the  transmission  loss  is 
shown  to  change  by  2  dB  as  the  thermocline  rises.  Calculations 
using  a  range-independent  broadband  parabolic  equation  model 
are  shown  to  produce  mean  transmission  loss  results  consistent 
with  experimental  observations  at  both  short  and  long  ranges. 
The  results  suggest  that  the  enduring  effects  of  passing  non¬ 
linear  internal  waves  on  mean  acoustic  transmission  loss  should 
be  observable,  significant,  and  predictable.  In  addition,  the  re¬ 
sults  have  implications  on  studies  such  as  geoacoustic  inversion 
using  transmission  loss  data. 

The  data  reported  in  this  paper  were  collected  during  the 
2006  Shallow  Water  Experiment  (SW06)  performed  on  the  con¬ 
tinental  shelf  off  the  coast  of  New  Jersey  [12].  Oceanographic 
results  are  described  in  Section  II.  Thermistor  data  show  how 
the  high-gradient  region  of  the  thermocline  is  depressed  by  non¬ 
linear  internal  waves  and  rises  only  gradually  to  its  background 
state.  Results  are  shown  to  be  consistent  with  an  oceanographic 
model  based  on  a  solution  to  the  Dubriel-Jacotin-Long  (DJL) 
equation.  Acoustic  results  are  described  in  Section  III.  Seven 
hours  of  acoustic  data  collected  in  the  1.5-10.5-kHz  band  at  a 
range  550  m  are  presented  for  a  period  before,  during,  and  im¬ 
mediately  after  the  passage  of  an  internal  wave.  Subsequent  to 
the  passage  of  the  wave,  another  eight  hours  of  data  were  col¬ 
lected  on  a  tow  track  with  the  range  steadily  increasing  from 
100  m  to  8  km.  Good  model/data  agreement  is  demonstrated 
for  both  fixed-  and  towed- source  scenarios.  The  results  are  dis¬ 
cussed  in  Section  IV. 

II.  Ocean  Data  and  Modeling 

Fig.  1  shows  the  effect  of  the  internal  tide  on  sound  speed  as 
observed  during  SW06.  The  sound  speed  is  plotted  as  a  function 
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Fig.  1.  Effect  of  the  slow  rising  of  the  thermocline  after  passing  nonlinear  in¬ 
ternal  waves.  Sound  speed  as  function  of  depth  for  a  48-h  period  beginning 
00:00:00  Universal  Time  Zone  (UTC)  August  18,  2006,  during  SW06.  Non¬ 
linear  internal  waves  depress  isotachs  at  approximately  the  M2  tidal  cycle.  Iso- 
tachs  gradually  rise  before  being  again  depressed  by  next  wave.  Superimposed 
are  regions  where  two  acoustic  data  sets  were  collected:  source  at  fixed  station 
(FS)  with  range  550  m;  and  towed  source  (TS).  (Mooring  data  were  acquired  by 
the  Woods  Hole  Oceanographic  Institution.) 


Fig.  2.  Positions  of  deployed  assets.  The  acoustic  receiving  array’s  position 
defines  the  origin.  Five  thermistor-chain  moorings  used  in  current  analysis  are 
labeled  14-17  and  54.  (Thermistor  chain  54  used  to  generate  Fig.  1.)  The  R/V 
Knorr  position  is  shown  for  fixed  station  (550  m)  acoustic  transmissions.  Po¬ 
sition  of  the  R/V  Oceanus  shown  at  21:14:00  UTC  August  19,  2006,  when 
nonlinear  internal  wave  labeled  NFIW  was  encountered.  Beginning  of  the  R/V 
Knorr' s  track  for  subsequent  towed  source  transmissions  is  also  shown. 


of  depth  for  48  h  [13].  For  this  section  of  data,  the  thermocline’ s 
high-gradient  region  was  thrust  several  meters  downward  as  in¬ 
ternal  waves  passed  at  approximately  the  M2  tidal  period  of  12.4 
h.  After  a  wave  has  passed,  the  thermocline  rises  only  gradually 
to  what  could  loosely  be  considered  a  background  profile.  Con¬ 
sequently,  while  the  main  part  of  the  internal  wave  may  be  near 
a  sensor  for  only  a  short  period,  the  complete  internal  tide  effect 
may  last  for  several  hours. 

The  pattern  displayed  in  Fig.  1  was  observed  throughout 
SW06.  The  experiment  featured  a  total  of  62  acoustics  and 
physical  oceanographic  moorings  deployed  in  a  geometry 
shaped  like  a  “T”  [12].  Of  present  interest  are  five  oceano¬ 
graphic  moorings  positioned  as  shown  in  Fig.  2  on  an  80-m 
isobath.  Additional  oceanographic  data  were  collected  from  the 
R/V  Knorr  and  R/V  Oceanus. 


IEEE  JOURNAL  OF  OCEANIC  ENGINEERING,  VOL.  35,  NO.  1,  JANUARY  2010 

During  the  48-h  period  shown  in  Fig.  1,  two  acoustic  data 
sets  were  collected.  In  the  first,  acoustic  signals  were  transmitted 
from  the  R/V  Knorr  at  a  fixed  station  to  a  vertical  receiving  array 
550  m  away.  In  the  second,  the  R/V  Knorr  slowly  towed  the 
acoustic  source  away  from  the  receiving  array.  Fig.  1  shows  the 
time  periods  when  these  two  acoustic  data  sets  were  collected 
and  Fig.  2  shows  the  positions  of  the  acoustic  assets  relative  to 
the  oceanographic  moorings.  In  this  section,  an  oceanographic 
model  is  developed  that  can  be  used  to  explain  the  acoustic  ob¬ 
servations  detailed  in  Section  III.  Note  that  the  oceanographic 
moorings  do  not  coincide  with  either  of  the  acoustic  tracks.  Con¬ 
sequently,  the  extent  to  which  these  relatively  distant  oceano¬ 
graphic  measurements  can  be  used  to  model  conditions  along 
the  acoustic  tracks  must  be  tested. 

The  most  dramatic  feature  of  Fig.  1  is  the  strong  depression 
in  the  thermocline  when  nonlinear  internal  waves  pass.  The  R/V 
Oceanus  was  positioned  as  shown  in  Fig.  2  at  21:14:00  UTC 
on  August  18,  2006,  when  it  was  passed  by  a  large  amplitude 
wave.  X-band  radar  measurements  from  the  ship  indicated  the 
wave’s  bearing  as  288°  and  speed  as  0.89  m/s  [14].  Using  these 
direct  measurements  as  ground  truth,  the  first  task  is  to  develop 
a  theoretical  model  for  the  rapid  thermocline  depression  caused 
by  the  internal  wave.  The  second  task  is  to  develop  a  practical 
model  that  uses  as  input  only  data  collected  on  the  oceano¬ 
graphic  moorings. 

A.  DJL  Model  for  Nonlinear  Internal  Waves 

Fully  nonlinear  solitary  waves  in  a  continuously  stratified 
flow  can  be  described  by  the  DJL  equation  [2],  [15].  A  soli¬ 
tary  internal  wave  is  a  single  wave  that  propagates  with  an 
unchanging  shape.  Unlike  other  descriptions  of  these  waves, 
no  small-amplitude  assumption  is  made;  quite  moderate  am¬ 
plitudes,  such  as  those  encountered  in  SW06,  exceed  the  small 
amplitude  requirements  of  these  other  approaches. 

Although  the  waves  of  present  concern  are  parts  of  a  wave 
train,  it  has  been  found,  both  in  SW06  and  in  other  measure¬ 
ments  [16],  that  the  DJL  equation  gives  a  rather  good  descrip¬ 
tion  of  the  first  wave  in  a  wave  train.  This  is  because  the  first 
wave  is  well  separated  from  the  later  ones  and  therefore,  it  can 
be  treated  as  a  solitary  wave.  For  a  given  background  stratifica¬ 
tion  that  exists  before  a  wave  arrives,  there  are  one-parameter 
families  of  solitary  waves.  Each  family  has  a  different  modal 
character.  In  the  Turkington  et  al.  [17]  iterative  algorithm  for 
finding  the  lowest  mode  family  of  solutions  of  the  DJL  equa¬ 
tion,  the  free  parameter  is  chosen  to  be  the  potential  energy  of 
the  solution.  Given  the  potential  energy,  all  other  characteris¬ 
tics  of  the  wave  are  determined.  To  find  which  DJL  solution  ap¬ 
plies  to  a  particular  wave,  some  characteristic,  normally  the  am¬ 
plitude,  is  fitted  to  the  wave.  From  the  corresponding  solution, 
other  variables  such  as  the  wave  speed  are  determined.  Present 
calculations  use  a  slight  modification  of  the  Turkington  et  al.  al¬ 
gorithm;  the  modification  improves  the  convergence  when  mea¬ 
sured  stratification  is  used. 

The  background  stratification  used  as  input  to  the  DJL  equa¬ 
tion  was  from  a  conductivity-temperature-depth  (CTD)  cast 
taken  on  August  1 8, 2006,  at  17:01 :00  UTC  from  the  R/V  Knorr, 
fully  4  h  before  the  next  major  internal  wave  was  encountered. 


YANG  et  al. :  EFFECT  OF  THE  INTERNAL  TIDE  ON  ACOUSTIC  TRANSMISSION  LOSS  AT  MIDFREQUENCIES 


5 


Fig.  3.  Temperature  record  at  depth  25  m  from  the  five  thermistor-chain  moor¬ 
ings  shown  in  Fig.  2.  (a)  Complete  48-h  period,  (b)  Enlargement  with  arrows 
pointing  to  arrival  of  the  large  nonlinear  internal  wave  for  each  sensor  [boxed 
region  in  (a)].  Observations  used  to  estimate  speed  and  bearing  of  wave. 


The  wave  speed  of  the  internal  wave  is  calculated  as  a  func¬ 
tion  of  the  internal  wave  amplitude.  Taking  the  wave  amplitude 
from  measurements  at  mooring  54  (Fig.  2),  the  predicted  wave 
speed  is  0.86  m/s.  The  conclusion  is  that  the  wave  speed  can  be 
predicted  from  first  principles  given  knowledge  of  the  ambient 
stratification  and  a  reasonable  estimate  of  the  wave  amplitude. 

B.  Mooring-Based  Model 

The  five  oceanographic  moorings  in  Fig.  2  sampled  the  tem¬ 
perature  every  30  s.  Moorings  14-17  had  three  sensors  each 
while  mooring  54  had  10.  Fig.  3(a)  shows  the  temperature  record 
from  each  mooring  at  approximate  depth  25  m  for  the  com¬ 
plete  48-h  record.  Fig.  3(b)  is  an  enlargement  emphasizing  the 
time  period  near  when  the  internal  wave  observed  by  the  RJV 
Oceanus  arrived  at  the  arrays.  The  arrows  indicate  when  the 
peak  internal  wave  displacement  was  observed  on  each  sensor. 

To  estimate  properties  of  the  internal  wave  using  only  the 
mooring  data,  the  wave  is  modeled  locally  as  a  plane  wave. 
Using  the  arrival  times  from  Fig.  3(b)  as  input,  the  bearing  and 
speed  of  the  wave  are  then  estimated  by  a  least  squares  fit.  The 
resulting  calculation  gives  speed  0.81  m/s  and  bearing  292° .  The 
agreement  with  ground- truth  speed  0.89  m/s  and  bearing  288° 
is  good.  The  moorings  were  in  water  depth  slightly  shallower 
than  the  RJV  Oceanus  so  a  slightly  lower  wave  speed  would  be 
expected. 

The  plane- wave  model  for  nonlinear  internal  waves  is  clearly 
limited.  SW06  was  performed  in  a  region  where  new  internal 
waves  were  spawned  rapidly  [18].  Wavefront  curvature  and 
merging  between  different  waves  may  be  significant.  Still,  the 
agreement  between  the  mooring-based  model  and  ground-truth 
measurements  is  encouraging;  it  suggests  that  the  mooring 
data  can  be  used  to  infer  at  least  the  gross  characteristics  of 
the  internal  wave  field  as  would  have  been  encountered  on  the 
acoustic  tracks  several  kilometers  away  (Fig.  2). 

Using  the  bearing  and  speed  estimated  from  the  plane- wave 
ocean  model,  it  can  be  determined  that  the  internal  wave  takes 
approximately  16  min  to  propagate  from  the  acoustic  receiving 
array  to  mooring  54.  Similarly,  it  takes  10  min  for  the  internal 


wave  to  propagate  from  the  RJV  Knorr  to  the  receiving  array 
for  the  550-m  range  transmission  experiment.  Based  on  these 
calculations,  the  following  ocean  model  is  proposed  for  use 
in  the  acoustic  simulations:  mooring  54  temperature  data  are 
smoothed  over  a  sliding  window  a  minimum  of  10  min  in  du¬ 
ration.  Then,  the  smoothed  temperature  records  are  offset  by  a 
minimum  of  16  min  to  correct  for  the  propagation  time  from 
the  acoustic  receiving  array.  Then,  the  smoothed,  offset  temper¬ 
ature  profiles  are  combined  with  salinity  data  to  calculate  the 
sound-speed  profile  (SSP)  at  the  acoustic  receiving  array.  This 
SSP  is  then  used  in  range-independent  acoustic  simulations  as 
detailed  in  the  following  section.  The  proposed  slowly  varying, 
range-independent  ocean  model  would  be  expected  to  fail  in  the 
immediate  vicinity  of  nonlinear  internal  waves.  However,  calcu¬ 
lations  in  the  following  section  will  show  it  to  be  adequate  for 
capturing  the  gross  effects  of  the  rising  thermocline  as  observed 
in  the  acoustic  experiment. 

III.  Acoustic  Data  and  Modeling 

As  noted  in  Section  II,  the  present  analysis  concerns  two 
acoustic  data  sets  collected  on  August  18-19, 2006.  For  the  first, 
the  RJV  Knorr  was  at  a  fixed  station,  550  m  from  the  receiving 
array.  The  550-m  range  was  selected  because  it  was  expected  to 
be  comparable  to  the  typical  width  of  a  nonlinear  internal  wave. 
At  this  range,  it  was  expected  that  the  different  acoustic  arrivals 
could  be  separated  from  one  another  and  studied  individually. 
For  the  second,  the  acoustic  source  was  slowly  towed  out  to  a 
maximum  range  of  8.1  km.  The  two  data  sets  allow  the  effects 
of  the  rising  thermocline  on  acoustic  propagation  to  be  exam¬ 
ined  over  a  range  of  temporal  and  spatial  scales. 

For  both  data  sets,  acoustic  signals  were  recorded  on  a 
moored  receiving  system  positioned  at  39°  01. 47'  N,  73° 
02.262'  W  (Fig.  2).  The  system  [19]  included  two  vertical 
subarrays,  each  with  four  elements:  a  shallow  subarray  with 
elements  at  depths  25.0,  25.2,  25.5,  and  26.4  m,  and  a  deep 
subarray  with  elements  at  50.0,  50.2,  50.5,  and  51.4  m.  Signals 
were  transmitted  using  an  ITC-2015  transducer  (International 
Transducer  Corporation)  positioned  at  nominal  depth  40  m  off 
the  stern  of  the  RJV  Knorr.  Linear  frequency-modulated  (LFM) 
chirp  signals,  20  ms  in  duration,  were  transmitted  with  an 
approximate  repetition  rate  of  19  s.  The  midfrequency  chirps 
swept  from  1.5-10.5  kHz  with  a  raised  cosine  window  and 
10%  taper.  In  subsequent  processing,  the  received  signals  were 
pulse  compressed  using  replicas  obtained  from  calibration  data. 
Details  of  how  the  replicas  were  generated  are  given  in  the 
Appendix. 

A.  Fixed  Source 

On  August  18,  2006,  beginning  at  15:16:00  UTC  and  ending 
at  22:32:00  UTC,  acoustic  data  were  collected  at  a  nominal 
range  550  m.  From  Fig.  1,  this  coincides  with  a  period  before, 
during,  and  immediately  after  the  passage  of  a  nonlinear  internal 
wave.  The  current  analysis  emphasizes  the  period  before  and 
after  the  wave;  the  period  during  the  wave  was  studied  previ¬ 
ously  [14]. 

A  total  of  1400  acoustic  transmissions  were  recorded  at  range 
550  m.  Fig.  4  shows  the  signals  received  at  depths  25  and  50  m 
after  pulse  compression  for  the  duration  of  the  experiment.  The 
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Fig.  4.  Acoustic  arrival  structure  with  source  at  fixed  station,  range  550  m.  Pulse  compression  output  versus  geotime,  (a)  Receiver  at  depth  25  m.  (b)  Receiver  at 
depth  50  m  (35-dB  dynamic  range).  Superimposed  are  simulation  results;  see  text  for  details. 


signals  were  aligned  using  the  bottom-bounce  arrival  as  refer¬ 
ence;  reduced  time  in  the  figure  is  with  reference  to  this  arrival. 
The  bottom  bounce  was  observed  to  be  the  most  stable  arrival 
throughout  the  entire  data  set.  The  figure  shows  the  first  30  ms 
of  the  compressed  signal,  sufficient  to  show  the  first  four  arrival 
groups  and  all  acoustic  paths  with  at  most  one  reflection  off  the 
sea  surface  and  one  reflection  off  the  seabed. 

There  are  features  in  Fig.  4  worth  considering  in  more  de¬ 
tail.  At  depth  25  m,  the  first  group  of  arrivals,  between  — 10  and 
—5  ms,  consists  of  three  paths:  two  waterborne  direct  paths  and  a 
sea  surface  reflection.  This  first  group  of  arrivals  gets  more  com¬ 
pact,  i.e.,  the  relative  arrival  time  difference  between  the  three 
paths  decreases  between  hours  15  and  19.  Shortly  after  hour  21, 
there  is  a  strong  acoustic  effect  produced  by  the  passage  of  an 
internal  wave.  It  takes  10  min  for  the  main  part  of  the  internal 
wave  to  pass  between  the  source  and  the  receiver  [14].  After  the 
internal  wave  passes,  the  arrival  time  and  the  time  width  of  the 
first  group  become  more  like  what  was  observed  around  hour 
15.  However,  gradual  changes  in  the  arrival  structure  between 
hours  15  and  19  are  less  apparent  than  at  depth  25  m. 

The  experimental  results  in  Fig.  4  demonstrate  that  the  ef¬ 
fects  produced  by  a  nonlinear  internal  wave  on  acoustic  signals 
are  manifest  for  several  hours.  The  primary  goal  of  the  acoustic 
modeling  is  to  reproduce  the  gradual  changes  in  the  acoustic 
arrival  structure  that  are  observed  in  the  data  after  the  wave 
has  passed.  Of  present  interest  is  the  slowly  varying  acoustic 
intensity  for  the  multiple  acoustic  arrivals.  Rapid  fluctuations 
in  intensity — scintillation — are  beyond  the  scope  of  this  study. 
Given  that  the  interest  is  only  in  the  slowly  varying  mean  in¬ 
tensity,  the  modeling  approach  is  to  assume  a  time- varying  but 
range-independent  representation  for  the  ocean  (Section  II).  A 
secondary  goal  is  to  see  if  this  relatively  simple  ocean  model 
can  account  for  the  acoustic  observations.  If  it  can,  it  implies 
that  bulk  acoustic  characteristics  can  be  predicted  successfully 
with  relatively  sparse  environmental  information  as  input. 

Acoustic  propagation  was  simulated  using  the  parabolic 
equation  (PE)  method  [20],  [21].  This  work  focuses  on 
data/model  comparison  in  the  frequency  band  of  1.5-6  kHz. 
Broadband  pulses  were  generated  by  Fourier  synthesis;  to  fill 
the  1.5-6-kHz  band,  1798  separate  single-frequency  PE  cal¬ 
culations  were  made.  The  sediment  parameters  (fluid  bottom) 


Fig.  5.  Simulated  acoustic  arrival  structure  out  to  range  550  m  for  receiver 
depth  25  m.  Beyond  range  450  m,  multiple  direct  and  surface  bounce  paths  are 
not  separated  in  time. 


used  were  as  follows:  sound  speed  1620  m/s  [22],  density 
1.85  g/cm3,  and  attenuation  0.5  dB/A  [23]. 

Fig.  5  shows  the  synthesized  time  series  for  receiver  depth 
25  m  as  a  function  of  range.  The  calculation  is  taken  to  the 
550-m  range  used  in  the  experiment  and  the  reduced  time  is 
again  referenced  to  the  bottom-bounce  path.  The  SSP  used  in  the 
calculation  is  derived  from  a  CTD  cast  taken  from  the  RJVKnorr 
at  17:01:00  UTC.  The  time  window  50  ms  is  longer  than  used 
in  Fig.  4  and  sufficient  to  include  two  additional  acoustic  paths: 
the  surface-bottom-surface  bounce  path  and  the  bottom-sur¬ 
face-bottom  bounce  path.  As  the  range  increases  beyond  about 
450  m,  the  early  arriving  paths  merge  and  form  a  group.  This 
first  group,  between  —10  and  —5  ms,  consists  of  three  arrivals. 
These  three  arrivals  correspond  to  a  fast  direct,  a  surface  bounce, 
and  a  slow  direct  arrival.  Beamforming  results  indicate  that,  in 
a  ray  picture,  all  three  arrivals  come  from  above  the  receiver. 

The  simulation  results,  using  the  CTD  cast  as  input  in  Fig.  5  at 
range  550  m,  can  be  compared  to  the  experimental  observations. 
The  simulation  results  are  overlaid  in  Fig.  4  by  the  solid  line 
(magenta  for  color  online)  at  hour  17.  For  both  depths  25  and 
50  m,  the  simulations  are  in  good  agreement  with  data  predicting 
both  number  of  arrivals  and  corresponding  arrival  times. 


YANG  et  al. :  EFFECT  OF  THE  INTERNAL  TIDE  ON  ACOUSTIC  TRANSMISSION  LOSS  AT  MIDFREQUENCIES 


7 


To  examine  the  thermocline  rising  effect  and  make  model/ 
data  comparisons  as  time  evolves,  it  is  necessary  to  use  time-de- 
pendent  SSPs  as  input  to  the  model.  However,  the  number  of 
CTD  casts  from  the  RJV  Knorr  was  limited  and  it  is  necessary 
to  use  as  input  data  collected  on  the  nearby  environmental  moor¬ 
ings.  Using  the  algorithm  developed  in  Section  II,  mooring  54 
data  are  used  to  construct  the  SSPs.  A  10-min  sliding  window 
was  applied  to  the  temperature  data.  The  measurements  were 
offset  by  20  min  to  compensate  for  the  travel  time  between  the 
midpoint  of  the  acoustic  track  and  mooring  54.  A  total  of  22  such 
profiles  were  generated  representing  the  changing  environment 
over  the  7-h  fixed  position  experiment. 

Simulation  results  at  selected  times  are  superimposed  over 
the  experimental  data  in  Fig.  4  as  well.  Note  that  the  simula¬ 
tion  curves  at  times  other  than  hour  17  (CTD  input)  use  the 
mooring  data  (white  lines).  Starting  once  more  with  depth  25  m, 
the  model/data  comparison  shows  good  agreement  not  only  in 
the  arrival  time  but  also  in  the  varying  compactness  of  the  first 
group  of  arrivals.  Similar  comparisons  of  the  later  arrivals  show 
similar  good  agreement.  For  depth  50  m,  the  agreement  between 
data  and  modeling  is  also  consistent. 

With  the  arrival  structure  well  characterized,  the  next  test  for 
the  model  is  signal  intensity.  As  noted  earlier,  the  first  group  of 
arrivals  for  depth  25  m  is  the  most  sensitive  to  the  rising  ther¬ 
mocline.  The  model  should  reproduce  the  observed  increase  in 
intensity  as  the  thermocline  rises.  It  should  also  show  a  drop  in 
intensity  as  the  thermocline  is  depressed,  the  most  dramatic  in¬ 
stance  being  during  the  passage  of  an  internal  wave.  This  can  be 
simply  explained  as  follows:  when  the  thermocline  is  depressed 
by  the  passing  wave,  the  receiver  at  25  m  can  be  regarded  as  out 
of  the  sound  channel;  as  the  thermocline  rises,  the  receiver  at 
25  m  is  then  in  the  channel,  and  therefore,  the  received  signal 
has  much  higher  intensity.  In  addition,  the  change  of  thermo¬ 
cline  depth  is  expected  to  affect  paths  that  go  through  the  mixed 
layer  more  than  the  other  paths  as  seen  in  the  first  group  of  ar¬ 
rivals  received  at  25  m. 

Fig.  6  shows  the  signal  intensity  (with  mean  removed)  of 
the  first  arrival  group  at  depth  25  m  for  the  entire  7-h  period. 
The  thin  line  shows  the  raw  data  and  the  heavy  line  is  the 
data  averaged  over  a  10-min  sliding  window.  Consistent  with 
Fig.  4,  the  correlation  between  increasing  signal  intensity  and 
the  rising  thermocline  is  apparent  between  hours  15  and  20. 
Between  hours  20  and  21,  the  thermocline  is  slightly  depressed 
followed  by  strong  depression  after  hour  21  when  the  internal 
wave  passes.  Signal  intensity  reaches  its  lowest  level  when  the 
thermocline  is  greatly  depressed  by  the  internal  wave.  Com¬ 
paring  measured  intensity  before  and  after  the  passage  of  the 
internal  wave  shows  approximately  5-dB  difference.  Simulation 
results  in  Fig.  6  are  shown  as  interconnected  dots.  The  acoustic 
simulations  reproduce  the  broad  characteristics  of  the  data  with 
particularly  good  agreement  before  hour  20.  Between  hour  20 
and  22,  the  simulation  over  predicts  the  observed  intensity.  This 
is  not  surprising  as  the  range-independent  environmental  model 
used  in  the  simulation  is  clearly  inadequate  when  the  nonlinear 
internal  waves  are  nearby.  The  internal  waves  drive  acoustic 
energy  into  the  seabed  and  increase  loss,  a  factor  not  included 
in  the  model.  At  hour  22,  after  the  internal  wave  has  passed,  the 
model/data  agreement  improves  as  might  be  expected. 


Fig.  6.  Model/data  comparison  for  integrated  signal  intensity  of  the  first  arrival 
group  at  depth  25  m.  Thin  line:  raw  data.  Thick  line:  lowpass  filtered  data.  Dotted 
line:  simulation. 
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Fig.  7.  Acoustic  arrival  structure  for  towed  source.  Pulse  compression  output 
versus  range  at  depth  25  m  using  frequency  band  1 .5-6  kHz  for  the  entire  tow 
track. 

B.  Towed  Source 

On  August  19,  2006,  beginning  at  00:30:00  UTC  and  ending 
at  09:03:00  UTC,  towed  source  acoustic  data  were  collected. 
The  tow  began  approximately  2  h  after  the  fixed  source  data 
presented  in  Section  III-A  and  3  h  after  the  passage  of  the  last 
internal  wave  event  (Fig.  1).  The  source  was  towed  by  the  RJV 
Knorr  at  depth  40  m  and  speed  0.26  m/s  along  an  80-m  isobath. 
The  range  to  the  acoustic  receiving  array  increased  steadily  from 
104  m  to  8.1  km;  Fig.  2  shows  the  beginning  of  the  track.  Using 
the  same  LFM  chirp  signals  as  used  earlier,  1285  transmissions 
were  recorded. 

Fig.  7  shows  the  complete  8 -km  towed- source  data  after  pulse 
compression  for  the  shallow  receiver  at  25  m.  In  this  figure,  the 
reduced  time  is  the  difference  between  actual  signal  propaga¬ 
tion  time  and  propagation  time  at  a  reference  sound  speed  of 
1495  m/s.  The  multiple  arrivals,  bouncing  between  the  surface 
and  the  bottom,  are  apparent  until  they  bundle  together  to  form 
distinct  groups  as  modes  at  approximately  4-km  range. 
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Fig.  8.  Reduced  transmission  loss  versus  range  for  towed  source,  (a)  Receiver  depth  25  m.  (b)  Receiver  depth  50  m. 


At  very  short  range,  the  first  three  arrivals  are  the  direct, 
surface,  and  bottom  bounce  paths.  As  range  increases  beyond 
0.4  km,  the  direct  and  surface  bounce  paths  merge  as  observed 
at  fixed  range  (Fig.  4)  and  predicted  in  numerical  simulation 
(Fig.  5).  Another  interesting  feature  is  the  short  paired  arrivals 
observed  repeatedly  for  range  up  to  4  km  between  —10  and 
20  ms.  They  appear  regularly  at  certain  ranges  and  last  about 
10  ms.  Using  ray  tracing  technique,  these  arrivals  are  identified 
as  ray  paths  that  are  reflected  off  the  bottom  and  turn  at  the  end 
of  the  warm  mixed  surface  layer  around  20  m.  One  of  the  pair 
of  arrivals  corresponds  to  launch  angles  going  upward  from  the 
source  (40  m)  to  the  receiver  (25  m)  and  the  other  corresponds 
to  launch  angles  going  downward.  The  widths  of  these  short  ar¬ 
rivals,  in  terms  of  range  and  time,  are  determined  by  how  wide 
the  ray  bundle  spans  at  receiver  depth  25  m. 

The  reduced  transmission  loss  for  the  shallow  (25  m)  and 
deep  (50  m)  receivers  is  shown  in  Fig.  8.  Results  are  plotted  at 
5  frequencies,  each  with  1-kHz  bandwidth.  In  this  context,  the 
reduced  transmission  loss  corrects  for  the  cylindrical  spreading 
and  the  frequency-dependent  water  absorption.  The  calibration 
data  discussed  in  the  Appendix  are  used  to  normalize  the 
tow  data  at  each  frequency  as  well.  Results  shown  in  Fig.  8 
have  been  smoothed  over  range  (125  m),  which  removes  the 
large  frequency-dependent  multipath  interference.  Without 
the  smoothing,  results  would  differ  greatly  between  different 
frequencies. 

Several  observations  can  be  made  from  Fig.  8.  For  both 
depths,  the  reduced  transmission  loss  generally  increases  with 
range.  The  shallow  receiver  exhibits  4-5  dB  more  loss  than  the 
deep  at  8  km.  The  extra  4-5-dB  loss  for  the  shallow  receiver,  in 
terms  of  normal  modes,  is  due  to  the  more  attenuative  modes 
residing  at  depth  25  m.  For  the  shallow  receiver,  the  reduced 
transmission  loss  is  only  weakly  dependent  on  frequency.  The 
loss  actually  decreases  by  approximately  2  dB  between  5 
and  7  km  before  again  increasing  with  range.  For  the  deep 
receiver,  the  reduced  transmission  loss  has  two  “plateaus”  at 
the  beginning  and  the  end  with  a  transition  region  in  the 
middle  between  4  and  6  km.  The  results  at  five  selected 
frequencies  are  fairly  close  in  the  two  plateau  regions  while 
there  is  a  slight  increase  in  loss  between  5  and  6  km  at  high 
frequencies,  e.g.,  at  4.5  kHz. 


Fig.  9  shows  a  sampling  of  SSPs  measured  at  mooring  54 
(Fig.  2)  during  the  towed  source  experiment.  Each  profile  repre¬ 
sents  an  average  over  a  time  window  20  min  in  duration.  Noted 
next  to  each  profile  is  the  time  at  the  center  of  the  window  and 
the  associated  range  between  the  towed  source  and  the  acoustic 
receiving  array.  The  double  sound-speed  duct,  discussed  rela¬ 
tive  to  Fig.  7(b),  is  apparent  within  the  first  3  h  and  diminishes 
as  the  thermocline  rises.  As  manifestations  of  the  rising  thermo- 
cline,  the  deep  sound  channel  axis  shifts  from  depth  45  to  37  m 
while  the  mixed  surface  layer  depth  shrinks  from  approximately 
20  to  10  m  over  the  8-h  period. 

The  goal  of  the  modeling  effort  is  to  develop  a  model  that  cap¬ 
tures  the  essential  features  of  the  experimental  results  in  Figs.  7 
and  8.  Certain  gross  features,  such  as  the  reduced  transmission 
loss  generally  increasing  with  range,  can  be  explained  by  inter¬ 
action  with  the  lossy  seabed.  Other  more  detailed  features,  such 
as  the  split  arrivals  in  Fig.  7  and  the  transmission  loss  plateaus 
in  Fig.  8,  depend  on  the  detailed  features  of  the  SSP  in  the  water 
column.  From  Fig.  9,  it  is  apparent  that  the  towed  source  data 
were  collected  during  a  period  when  the  thermocline  was  rising. 
To  capture  the  detailed  feature  of  the  data,  the  model  must  in¬ 
clude  a  rising  thermocline. 

As  in  Section  III-A,  the  acoustic  model  assumes  a  range-in¬ 
dependent  ocean  with  the  slowly  varying  SSP  derived  from 
mooring  data.  Clearly,  this  is  an  approximation  that  becomes 
increasingly  coarse  as  the  source-receiver  range  increases  and 
the  towed  source  gets  further  from  the  oceanographic  moorings 
(Fig.  2).  A  total  of  27  SSPs  from  mooring  54,  each  a  20-min 
average,  represented  the  evolving  environment  over  8.5  h.  A 
particular  profile  was  used  in  the  acoustic  simulation  only  for 
source-receiver  ranges  applicable  for  that  20-min  period.  The 
final  modeling  result  is  obtained  by  piecing  together  the  27 
segments  of  reduced  transmission  loss  at  their  corresponding 
ranges.  Bottom  parameters  were  unchanged  from  Section  III-A. 

Fig.  10  compares  the  model  to  data  for  the  reduced  transmis¬ 
sion  loss  at  receiver  depths  25  and  50  m.  The  comparison  is  at 
2.5  kHz  with  1-kHz  bandwidth.  The  thin  line  represents  data 
while  the  thick  line  represents  simulation  results.  In  general,  for 
both  receiver  depths,  there  is  good  model/data  agreement.  As 
might  be  expected,  the  agreement  is  particularly  good  for  short 
ranges — less  than  perhaps  4  km — where  the  model  reproduces 
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Fig.  9.  Time-evolving  SSP  during  the  8-h  towed  source  experiment.  X-axis:  sound  speed  for  each  profile  in  range  [1490  1530]  m/s. 


Fig.  10.  Model/data  comparison  for  reduced  transmission  loss  versus  range,  towed  source,  frequency  2.5  kHz.  Thin  line:  lowpass  filtered  data;  thick  line:  simu¬ 
lation.  (a)  Receiver  depth  25  m.  (b)  Receiver  depth  50  m. 


even  the  finer  features.  At  depth  25  m,  the  model  captures  the 
plateau  between  5  and  7  km  and  the  subsequent  increase  in  re¬ 
duced  transmission  loss  at  greater  ranges.  The  largest  error  is 
for  the  50-m  receiver  at  ranges  between  6  and  7  km  where  the 
model  underpredicts  the  reduced  transmission  loss.  It  has  been 
observed  through  numerical  experiments  that  the  deep  receiver 
can  sometimes  exhibit  significant  convergence-zone-like  oscil¬ 
lations  due  to  the  complexity  of  the  SSP.  The  occurrence  of  such 
oscillations  depends  on  environmental  information  such  as  rel¬ 
ative  depth  between  source  and  the  deep  sound  channel  axis 
and  depth  of  the  shallow  sound  channel  axis.  The  dip  for  the 
deep  channel  happens  to  be  in  the  valley  of  that  type  of  oscil¬ 
lation.  Furthermore,  the  4-5-dB  difference  in  loss  between  the 
shallow  and  deep  receivers  is  also  well  characterized  by  simu¬ 
lation  results. 

IV.  Summary 

Nonlinear  internal  waves  depress  the  high-gradient  region  of 
the  thermocline.  After  the  waves  have  passed,  the  thermocline 
rises  only  gradually  towards  its  prewave  level.  Rapid  depression 


of  the  thermocline  followed  by  gradual  rising,  both  manifesta¬ 
tions  of  the  internal  tide,  was  observed  repeatedly  during  SW06. 
This  work  emphasizes  the  effect  that  the  gradually  rising  ther¬ 
mocline  has  on  acoustic  propagation  in  the  midfrequency  band. 
The  effect  is  shown  to  be  significant  for  the  receiver  within  the 
thermocline  (25  m):  at  fixed  range  (550  m),  the  arrival  structure 
changes  and  there  is  a  5-dB  change  in  the  intensity  of  the  first 
arrival  group.  Similarly,  the  towed  source  data  shows  a  2-dB  in¬ 
crease  in  total  intensity  as  the  thermocline  rises. 

Using  nearby  mooring  data,  a  simple  plane- wave  model  for 
the  internal  wave  was  developed.  The  speed  and  bearing  for  the 
internal  wave  produced  by  this  simple  model  were  shown  to  be 
in  good  agreement  with  both  concurrent  radar  observations  and 
theoretical  calculations  based  on  the  DJL  equation.  The  bearing 
and  speed  estimates  were  used  to  calculate  the  time  offset  be¬ 
tween  the  environment  as  measured  at  the  moorings  and  what 
would  have  been  encountered  along  the  acoustic  tracks.  A  model 
for  the  SSP  results  that  can  be  used  in  acoustic  simulations 
is  as  follows:  the  sound  speed  is  treated  as  range  independent 
but  slowly  varying  in  time  as  the  thermocline  rises.  Without  in¬ 
cluding  this  time  dependence,  the  5-dB  change  in  transmission 
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Fig.  11.  Calibration  data  for  120  consecutive  transmissions,  depth  25  m.  Raw 
LFM  signals  without  pulse  compression.  Superimposed  heavy  line  shows 
source-receiver  range  derived  from  gyroscope  and  GPS  data  and  used  to 
correct  relative  motion. 


loss  observed  in  Fig.  6  would  not  be  captured  by  the  model. 
The  time-dependent  model  also  captures  the  finer  features  in 
the  towed  source  data  (Fig.  10),  particularly  for  ranges  less  than 
4  km.  A  single  SSP  is  not  adequate  for  predicting  transmission 
loss  levels.  In  addition,  if  transmission  loss  data  is  used  for  geoa¬ 
coustic  inversion,  as  many  do,  it  is  important  to  consider  the  ef¬ 
fect  of  the  rising  thermocline. 

Internal  waves  are  one  manifestation  of  the  internal  tide.  An¬ 
other  manifestation  is  the  slowly  rising  thermocline  that  occurs 
after  the  wave  has  passed.  With  respect  to  acoustic  propagation, 
while  the  effect  of  internal  waves  may  be  more  dramatic,  the  ef¬ 
fect  of  the  rising  thermocline  may  be  much  longer  lived.  Results 
from  this  paper  show  the  latter  effect  to  be  both  observable  in 
the  data  and  predictable  with  a  simple  model. 

Appendix 

The  analysis  in  Section  III  uses  replica  to  time  compress  LFM 
signals.  Pulse  compression  has  been  widely  used  in  signal  pro¬ 
cessing  as  it  helps  to  achieve  the  desired  range  resolution  with  a 
reduced  power  of  the  transmitter.  This  Appendix  outlines  how 
the  replicas  are  generated.  Compensating  for  relative  motion  be¬ 
tween  the  acoustic  source  and  receiver  is  shown  to  be  an  impor¬ 
tant  step  in  generating  reliable  replicas. 

Calibration  data  were  taken  on  August  1 1 , 2006.  Signals  were 
transmitted  from  a  source  at  depth  30  m  off  the  stern  of  the 
RJV  Knorr  and  received  50  m  away  at  the  receiving  array.  Of 
present  interest  is  the  LFM  part  of  the  signal  sweeping  from  1.5 
to  10.5  kHz  over  20  ms.  Fig.  1 1  shows  the  LFM  part  of  the  signal 
received  at  depth  25  m  for  120  consecutive  transmissions  with 
repetition  rate  10  s.  Two  points  should  be  observed.  First,  the 
arrival  time  wanders  by  5  ms  over  the  duration  of  the  calibration. 
Second,  without  pulse  compression,  the  different  acoustic  paths 
are  not  separated  in  time.  The  direct  path  signal,  for  example, 
overlaps  with  the  surface  reflected  path. 

The  wander  in  Fig.  1 1  is  greater  than  what  can  be  attributed 
to  ocean  variability.  The  wander  instead  is  presumably  due  to 
relative  motion  between  the  source  deployed  off  the  ship  and 
the  moored  receiver.  Using  the  ship’s  global  positioning  system 
(GPS)  and  gyroscope  measurements  to  estimate  the  source-re¬ 


ceiver  range  independent  of  the  acoustic  observations  can  test 
this  presumption.  The  heavy  line  superimposed  on  Fig.  11 
shows  the  travel  time  calculated  using  the  estimated  source-re¬ 
ceiver  range  divided  by  the  mean  sound  speed.  The  agreement 
is  good  and  suggests  the  practicality  of  compensating  for 
source-receiver  motion  when  using  multiple  transmissions  to 
estimate  the  replica. 

Estimating  the  replica  signal  when  there  is  overlap  between 
multiple  acoustic  paths  involves  several  steps.  First,  the  signal 
is  pulse  compressed  to  separate  the  different  arrivals.  The  direct 
arrival  is  isolated  and  decompressed.  Then,  the  decompressed 
direct  arrivals  from  each  of  the  transmissions  are  aligned  after 
compensating  for  source-receiver  motion.  Coherent  averaging 
across  the  transmissions  yields  the  replica.  The  procedure  is  re¬ 
peated  for  each  element  in  the  receiving  array. 
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